1  Notre Mesure: Le Cultural Gender Normativity Index (CGNI)

code R
library(foreign)
library(questionr) 
library(ggplot2) 
library(tidyverse) 
library(ggmosaic)
library(GGally) 
library(dataMaid)
library(dplyr) 
library(GDAtools)
library(FactoMineR)
library(gtsummary) 
library(factoextra)
library(gtsummary) 
library(kableExtra)
library(RColorBrewer) 
library(FactoMineR) 
library(xtable) 
library(explor)
code R
data<-read.csv2("pc18_quetelet_octobre2023.csv")
code R
data$Sex <- factor(data$SEXE, 
                            levels = c(1, 2), 
                            labels = c("Men", "Women"))
my_data_frame <- data |> dplyr::rename( 
  Knitting = A1001 , 
  Cards_games = A1002,  
  Gambling = A1003 , 
  Cooking = A1004 , 
  DIY = A1005  ,
  Vegetable_garden = A1006 , 
  Ornamental_garden = A1007,  
  Fishing_hunting = A1008 , 
  Collection = A1009  ,
  Vehicle_custom = A1010 , 
  No_Amateur=A1011,
  Making_music = A1901  ,
  Diary = A1902  ,
  Writing = A1903,  
  Painting = A1904,  
  Montage = A1905 , 
  Circus = A1906  ,
  Pottery = A1907 , 
  Theater = A1908 , 
  Drawing = A1909 , 
  Dancing = A1910,  
  Photography = A1911  ,
  Genealogy = A1912 , 
  Science = A1913  ,
  None = A1914  ,
  Video_games = B1  ,
  TV = C1  ,
  Radio = E1  ,
  Library = F1  ,
  Museums = H112,  
  Internet = I4 , 
  Concert = G2413 )


my_data_frame$Video_games <- ifelse(my_data_frame$Video_games == 1, 1, 0)
my_data_frame$TV <- ifelse(my_data_frame$TV == 5, 0, 1)
my_data_frame$Radio <- ifelse(my_data_frame$Radio == 5, 0, 1)
my_data_frame$Library<- ifelse(my_data_frame$Radio == 1, 1, 0)
my_data_frame$Museums<- ifelse(my_data_frame$Museums == 1, 0, 1)
my_data_frame$Internet<- ifelse(my_data_frame$Internet == 5, 0, 1)
my_data_frame$Concert<- ifelse(my_data_frame$Concert == 1, 0, 1)


my_data_frame <- my_data_frame %>%
  mutate(satisfaction = case_when(
    A2 %in% 1 ~ "Low",      # 1 à 4 -> Low
    A2 %in% 2 ~ "Medium",   # 5 à 7 -> Medium
    A2 %in% 3 ~ "High",    # 8 à 10 -> High
    TRUE ~ NA_character_              
  ))  
my_data_frame <- my_data_frame %>%
  mutate(Income = case_when(
    CRITREVENU %in% 1:4 ~ "Low",      # 1 à 4 -> Low
    CRITREVENU %in% 5:7 ~ "Medium",   # 5 à 7 -> Medium
    CRITREVENU %in% 8:10 ~ "High",    # 8 à 10 -> High
    TRUE ~ NA_character_              
  ))

my_data_frame <- my_data_frame %>%
  mutate(Health = case_when(
    A15 %in% 1:2 ~ "Good",      # 1 à 4 -> Low
    A15 %in% 3 ~ "Medium",   # 5 à 7 -> Medium
    A15 %in% 4:5 ~ "Bad",    # 8 à 10 -> High
    TRUE ~ NA_character_              
  ))


my_data_frame <- my_data_frame %>%
  mutate(Couple = case_when(
    VITENCOUPLE %in% 1:2 ~ "Yes",     
    VITENCOUPLE %in% 3~ "No",  
    
    TRUE ~ NA_character_              
  ))

quartiles <- quantile(data$AGE, probs = c(0.25, 0.5, 0.75), na.rm = TRUE)

my_data_frame$age_group <- cut(
  my_data_frame$AGE,
  breaks = 4,  # Automatically divide into 4 slices
  labels = c("[15-38[", "[38-54[", "[54-67[", "[67-97["),  # Labels optionnels
  include.lowest = TRUE 
)
code R
table <- my_data_frame |>
  tbl_summary(
    include = c(  "Knitting" , 
                  "Cards_games",  
                  "Gambling" , 
                  "Cooking" , 
                  "DIY"  ,
                  "Vegetable_garden" , 
                  "Ornamental_garden",  
                  "Fishing_hunting" , 
                  "Collection"  ,
                  "Vehicle_custom", 
                  "Making_music"   ,
                  "Diary" ,
                  "Writing" ,  
                  "Painting",  
                  "Montage"  , 
                  "Circus"   ,
                  "Pottery" , 
                  "Theater" , 
                  "Drawing" , 
                  "Dancing",  
                  "Photography"  ,
                  "Genealogy" , 
                  "Science"  ,
                  "None"  ,
                  "No_Amateur",
                  "Video_games"  ,
                  "TV" ,
                  "Radio"  ,
                  "Library"  ,
                  "Museums",  
                  "Internet", 
                  "Concert"),
    by = "Sex"
  ) |>
  add_overall(last = TRUE) |>
  add_p()

table
Characteristic Men
N = 4,1621
Women
N = 5,0721
Overall
N = 9,2341
p-value2
Knitting 73 (1.8%) 1,370 (27%) 1,443 (16%) <0.001
Cards_games 1,899 (46%) 2,764 (54%) 4,663 (50%) <0.001
Gambling 990 (24%) 970 (19%) 1,960 (21%) <0.001
Cooking 1,685 (40%) 3,473 (68%) 5,158 (56%) <0.001
DIY 2,667 (64%) 2,283 (45%) 4,950 (54%) <0.001
Vegetable_garden 1,335 (32%) 1,245 (25%) 2,580 (28%) <0.001
Ornamental_garden 1,793 (43%) 2,198 (43%) 3,991 (43%) 0.8
Fishing_hunting 722 (17%) 212 (4.2%) 934 (10%) <0.001
Collection 395 (9.5%) 281 (5.5%) 676 (7.3%) <0.001
Vehicle_custom 264 (6.3%) 60 (1.2%) 324 (3.5%) <0.001
Making_music 1,312 (32%) 1,830 (36%) 3,142 (34%) <0.001
Diary 285 (6.8%) 1,206 (24%) 1,491 (16%) <0.001
Writing 410 (9.9%) 746 (15%) 1,156 (13%) <0.001
Painting 667 (16%) 1,303 (26%) 1,970 (21%) <0.001
Montage 843 (20%) 586 (12%) 1,429 (15%) <0.001
Circus 116 (2.8%) 179 (3.5%) 295 (3.2%) 0.044
Pottery 264 (6.3%) 705 (14%) 969 (10%) <0.001
Theater 483 (12%) 800 (16%) 1,283 (14%) <0.001
Drawing 864 (21%) 1,285 (25%) 2,149 (23%) <0.001
Dancing 410 (9.9%) 1,782 (35%) 2,192 (24%) <0.001
Photography 1,175 (28%) 1,176 (23%) 2,351 (25%) <0.001
Genealogy 513 (12%) 575 (11%) 1,088 (12%) 0.14
Science 611 (15%) 469 (9.2%) 1,080 (12%) <0.001
None 1,394 (33%) 1,261 (25%) 2,655 (29%) <0.001
No_Amateur 276 (6.6%) 359 (7.1%) 635 (6.9%) 0.4
Video_games 1,760 (42%) 1,827 (36%) 3,587 (39%) <0.001
TV 3,878 (93%) 4,806 (95%) 8,684 (94%) 0.001
Radio 3,522 (85%) 4,186 (83%) 7,708 (83%) 0.007
Library 3,522 (85%) 4,186 (83%) 7,708 (83%) 0.007
Museums 4,101 (99%) 5,009 (99%) 9,110 (99%) 0.4
Internet 3,459 (83%) 4,185 (83%) 7,644 (83%) 0.4
Concert 3,344 (80%) 4,234 (83%) 7,578 (82%) <0.001
1 n (%)
2 Pearson’s Chi-squared test
Encadré Technique 1: Détails de la méthode ACM

L’ACM permet d’explorer les relations entre plusieurs variables qualitatives en projetant les individus et les modalités dans un espace de faible dimension. Elle est souvent utilisée pour analyser des questionnaires et des tableaux de contingence complexes.

Principaux résultats d’une ACM

  1. Inertie totale

Mesure la dispersion des données et est donnée par :

\(I_{\text{total}} = \frac{q}{q-1} \sum_{k} \lambda_k\)

\(\lambda_k\) sont les valeurs propres et \(q\) est le nombre total de modalités.

  1. Valeurs propres \(\lambda_k\)

Elles indiquent la variance expliquée par chaque axe factoriel. Plus une valeur propre est élevée, plus l’axe correspondant est important dans l’analyse.

  1. Rapports de corrélation \(\eta^2\)

Le rapport de corrélation \(\eta^2\) mesure la liaison entre une variable et un axe factoriel :

\(\eta^2 = \frac{\sum_{i} f_i d_{i,k}^2}{\sum_{i} f_i d_{i}^2}\)

où $f_i $ est la fréquence de l’individu/modalité $i $, et \(d_{i,k}\)est sa distance à l’axe \(k\) .

  1. Coordonnées des individus et modalités

Elles sont obtenues à partir des vecteurs propres et permettent la représentation graphique des données :

\(C_{i,k} = \frac{v_{i,k}}{\sqrt{\lambda_k}}\)

\(v_{i,k}\) est le vecteur propre associé à l’axe \(k\).

  1. Cos² (Qualité de représentation)

Indique dans quelle mesure un point est bien représenté sur un axe donné. Une valeur proche de **1** signifie que la projection sur cet axe est pertinente.

  1. Contributions

Elles mesurent l’importance d’une modalité ou d’un individu dans la construction d’un axe. Plus une contribution est élevée, plus l’élément joue un rôle important dans l’interprétation de l’axe.

code R
pratiques_cols_1 <- c( "Knitting" , 
                     "Cards_games",  
                     "Gambling" , 
                     "Cooking" , 
                     "DIY"  ,
                     "Vegetable_garden" , 
                     "Fishing_hunting" , 
                     "Collection"  ,
                     "Vehicle_custom",
                     "Making_music"   ,
                     "Diary" ,
                     "Writing" ,  
                     "Painting",  
                     "Montage"  , 
                     "Pottery" , 
                     "Theater" , 
                     "Drawing" , 
                     "Dancing",  
                     "Photography"  ,
                     "Genealogy" , 
                     "Science"  ,
                     "None" ,
                     "Video_games"  ,
                     "Library"  ,
                     "Concert"
                            )


#MCA

# Add Sex Column to the selection
cols_of_interest_1 <- c("Sex", "AGE", pratiques_cols_1)

# Build a new dataframe with these columns
data_pratiques <- my_data_frame[, cols_of_interest_1]
data_pratiques$AGE <- cut(data_pratiques$AGE, 
                          breaks = quantile(data_pratiques$AGE, probs = seq(0, 1, 0.25), na.rm = TRUE), 
                          include.lowest = TRUE)
ra_data <- na.omit(data_pratiques)

cols_to_factor <- c(  "Knitting" , 
                      "Cards_games",  
                      "Gambling" , 
                      "Cooking" , 
                      "DIY"  ,
                      "Vegetable_garden" , 
                      "Fishing_hunting" , 
                      "Collection"  ,
                      "Vehicle_custom",
                      "Making_music"   ,
                      "Diary" ,
                      "Writing" ,  
                      "Painting",  
                      "Montage"  , 
                      "Pottery" , 
                      "Theater" , 
                      "Drawing" , 
                      "Dancing",  
                      "Photography"  ,
                      "Genealogy" , 
                      "Science"  ,
                      "None"  ,
                      "Video_games"  ,
                      "Library"  ,
                      "Concert")

# apply as.factor to these columnns
ra_data[cols_to_factor] <- lapply(ra_data[cols_to_factor], as.factor)

# running MCA with FactoMiner
acm2_fm <- ra_data |> 
  FactoMineR::MCA(
    ncp = Inf,
    graph = TRUE,
    quali.sup = 1:2
  )

code R
# Extract modality names
modalites_names <- rownames(acm2_fm$var$coord)

# Check modality names
head(modalites_names)
[1] "Knitting_0"    "Knitting_1"    "Cards_games_0" "Cards_games_1"
[5] "Gambling_0"    "Gambling_1"   
code R
# Extract coordinates for dimension 2
coord_dim2_modalites <- acm2_fm$var$coord[, 2]

# Create a table associating the modalities and their coordinates in dimension 2
modalites_coord <- data.frame(Modalite = modalites_names, Coord_Dim2 = coord_dim2_modalites)



# Keep only the two necessary columns
modalites_coord_selected <- modalites_coord[, c("Modalite", "Coord_Dim2")]

print(modalites_coord_selected)
                             Modalite  Coord_Dim2
Knitting_0                 Knitting_0  0.04314603
Knitting_1                 Knitting_1 -0.23295268
Cards_games_0           Cards_games_0 -0.21072530
Cards_games_1           Cards_games_1  0.20656774
Gambling_0                 Gambling_0 -0.17433328
Gambling_1                 Gambling_1  0.64698995
Cooking_0                   Cooking_0 -0.06432635
Cooking_1                   Cooking_1  0.05083253
DIY_0                           DIY_0 -0.59097111
DIY_1                           DIY_1  0.51145863
Vegetable_garden_0 Vegetable_garden_0 -0.29999865
Vegetable_garden_1 Vegetable_garden_1  0.77371744
Fishing_hunting_0   Fishing_hunting_0 -0.16671349
Fishing_hunting_1   Fishing_hunting_1  1.48150102
Collection_0             Collection_0 -0.05539977
Collection_1             Collection_1  0.70134797
Vehicle_custom_0     Vehicle_custom_0 -0.05128080
Vehicle_custom_1     Vehicle_custom_1  1.41022187
Making_music_0         Making_music_0  0.09403755
Making_music_1         Making_music_1 -0.18232870
Diary_0                       Diary_0  0.09378209
Diary_1                       Diary_1 -0.48702531
Writing_0                   Writing_0  0.06770001
Writing_1                   Writing_1 -0.47308017
Painting_0                 Painting_0  0.04361213
Painting_1                 Painting_1 -0.16081142
Montage_0                   Montage_0 -0.07125251
Montage_1                   Montage_1  0.38917132
Pottery_0                   Pottery_0  0.01201096
Pottery_1                   Pottery_1 -0.10244644
Theater_0                   Theater_0  0.07187594
Theater_1                   Theater_1 -0.44542916
Drawing_0                   Drawing_0  0.04570585
Drawing_1                   Drawing_1 -0.15068681
Dancing_0                   Dancing_0  0.14396501
Dancing_1                   Dancing_1 -0.46250072
Photography_0           Photography_0 -0.07539973
Photography_1           Photography_1  0.22074706
Genealogy_0               Genealogy_0 -0.01967874
Genealogy_1               Genealogy_1  0.14733732
Science_0                   Science_0 -0.04008198
Science_1                   Science_1  0.30261895
None_0                         None_0 -0.06572453
None_1                         None_1  0.16286315
Video_games_0           Video_games_0 -0.22185919
Video_games_1           Video_games_1  0.34927205
Library_0                   Library_0 -0.65936864
Library_1                   Library_1  0.13053925
Concert_0                   Concert_0 -0.23244943
Concert_1                   Concert_1  0.05079655
code R
# Initialize a vector to store the index of each individual
ra_data$indice_culturel <- 0

# Browse each individual
for (i in 1:nrow(ra_data)) {
  
  # Initialize individual's index to 0
  indice_individu <- 0
  
  # Browse each practice column (columns 3 to 27)
  for (pratique in 3:27) {
    
    # Retrieve the individual's response for this practice (0 or 1)
    reponse <- ra_data[i, pratique]
    
    # If the answer is 1, add the coordinate of the corresponding modality to the index.
    if (reponse == 1) {
      
      # Create the modality name (e.g. “knitting_1” or “knitting_0”)
      nom_modalite_1 <- paste0(names(ra_data)[pratique], "_1")
      nom_modalite_0 <- paste0(names(ra_data)[pratique], "_0")
      
      # Find the coordinate associated with the corresponding modality
      if (nom_modalite_1 %in% modalites_coord$Modalite) {
        indice_individu <- indice_individu + modalites_coord$Coord_Dim2[modalites_coord$Modalite == nom_modalite_1]
      }
      if (nom_modalite_0 %in% modalites_coord$Modalite) {
        indice_individu <- indice_individu + modalites_coord$Coord_Dim2[modalites_coord$Modalite == nom_modalite_0]
      }
    }
  }
  
  # Assign the calculated index to the individual
  ra_data$indice_culturel[i] <- indice_individu
}
####Normalisation

# Calculate minimum and maximum index values
min_indice <- min(ra_data$indice_culturel, na.rm = TRUE)
max_indice <- max(ra_data$indice_culturel, na.rm = TRUE)

# Normalize index
ra_data$indice_culturel_normalise <- (ra_data$indice_culturel - min_indice) / (max_indice - min_indice)

# Check results

head(ra_data[, c("indice_culturel", "indice_culturel_normalise")])
  indice_culturel indice_culturel_normalise
1       2.8024423                 0.8624142
2       0.5427446                 0.5073268
3      -0.4557477                 0.3504244
4      -0.3429253                 0.3681533
5       0.0409659                 0.4284777
6      -0.2367887                 0.3848315

Notre indice est donc construit de la façon suivante:

\[I_{1j} = \sum_{k=1}^{Z} w_{1k} \cdot X_{k j}\]

Dans cette expression, \(I_{1j}\) désigne l’indice de l’individu \(j\), tandis que \(w_{1k}\) représente le poids associé à chaque variable culturelle \(X_{kj}\). La somme englobe toutes les variables culturelles \(Z\), ce qui nous permet de saisir l’engagement culturel global de l’individu.

code R
my_data_frame$identity<-ra_data$indice_culturel_normalise
my_data_frame$indice<-ra_data$indice_culturel
ggplot(my_data_frame, aes(x = identity, fill = Sex)) +
  geom_density(alpha = 0.5) +
  scale_fill_manual(values = c("blue", "pink")) + 
  labs(title = "Density of The Normalized Cultural Index by Sexe",
       x = "Normalized Cultural Index",
       y = "Density",
       fill = "Sexe") +
  theme_minimal()

code R
stat_des_indice<-my_data_frame%>%
  tbl_summary(include=c("identity"),by = "Sex",
              statistic = list(
            all_continuous() ~ "{min} - {max}"
        )) %>%
  add_overall(last = TRUE) %>%
  add_p()
stat_des_indice
Characteristic Men
N = 4,1621
Women
N = 5,0721
Overall
N = 9,2341
p-value2
identity 0.01 - 1.00 0.00 - 0.91 0.00 - 1.00 <0.001
1 Min - Max
2 Wilcoxon rank sum test
code R
data$identity<- my_data_frame$identity
my_data_frame$identity<-data$identity
# List of cultural activities
cultural_activities <- c(
   "Knitting" , 
                      "Cards_games",  
                      "Gambling" , 
                      "Cooking" , 
                      "DIY"  ,
                      "Vegetable_garden" , 
                      "Fishing_hunting" , 
                      "Collection"  ,
                      "Vehicle_custom",
                      "Making_music"   ,
                      "Diary" ,
                      "Writing" ,  
                      "Painting",  
                      "Montage"  , 
                      "Pottery" , 
                      "Theater" , 
                      "Drawing" , 
                      "Dancing",  
                      "Photography"  ,
                      "Genealogy" , 
                      "Science"  ,
                      "None"  ,
                      "Video_games"  ,
                      "Library"  ,
                      "Concert"
)

# Create a data frame to store the results
result_table <- data.frame(Activity = character(0), Accuracy = numeric(0))

for (activity in cultural_activities) {
  # Perform a Probit regression for the current cultural activity
  model_formula <- as.formula(paste(activity, "~ identity"))
  model <- glm(model_formula, data = my_data_frame, family = binomial(link = "probit"))
 
  # Calculate predictions
  predicted <- ifelse(predict(model, type = "response") >= 0.5, 1, 0)
 
  # Calculate the accuracy
  correct_predictions <- sum(predicted == my_data_frame[[activity]])
  total_predictions <- nrow(my_data_frame)
  accuracy <- (correct_predictions / total_predictions) * 100
 
  # Add the result to the data frame
  result_table <- rbind(result_table, data.frame(Activity = activity, Accuracy = accuracy))
}
result_table
           Activity Accuracy
1          Knitting 84.37297
2       Cards_games 51.00715
3          Gambling 79.18562
4           Cooking 56.55187
5               DIY 51.25623
6  Vegetable_garden 72.89365
7   Fishing_hunting 93.95712
8        Collection 92.54927
9    Vehicle_custom 96.85943
10     Making_music 68.16114
11            Diary 85.85662
12          Writing 88.16331
13         Painting 78.76327
14          Montage 84.52458
15          Pottery 89.50617
16          Theater 86.53888
17          Drawing 76.73814
18          Dancing 79.66212
19      Photography 74.53974
20        Genealogy 88.21746
21          Science 88.30409
22             None 70.19710
23      Video_games 62.06411
24          Library 83.20338
25          Concert 81.92549
code R
my_data_frame$identity_scale <- cut(
  my_data_frame$identity,
  breaks = 7,  # Automatically divide into 7 slices
  labels = c("Very Feminine", "1", "2", "3", "4", "5", "Very Masculine"), 
  include.lowest = TRUE  
)
table1 <-
  my_data_frame |> 
  tbl_summary(include = c(identity_scale),
  by=Sex ,)|> 
    add_p()
table1
Characteristic Men
N = 4,1621
Women
N = 5,0721
p-value2
identity_scale

<0.001
    Very Feminine 9 (0.2%) 218 (4.3%)
    1 403 (9.7%) 1,506 (30%)
    2 2,139 (51%) 2,669 (53%)
    3 985 (24%) 576 (11%)
    4 509 (12%) 90 (1.8%)
    5 96 (2.3%) 10 (0.2%)
    Very Masculine 21 (0.5%) 3 (<0.1%)
1 n (%)
2 Pearson’s Chi-squared test

1.1 Indice Alternatif: Méthode de LASSO

code R
summary(my_data_frame$Library)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  1.0000  1.0000  0.8347  1.0000  1.0000 
code R
my_data_frame$SEXE<-as.factor(my_data_frame$SEXE)
reg_lm<- lm(Sex~ Knitting+ 
                  Cards_games+  
                  Gambling+ 
                  Cooking+ 
                  DIY+
                  Vegetable_garden+ 
                  Ornamental_garden+  
                  Fishing_hunting+ 
                  Collection+
                  Vehicle_custom + 
                  Making_music +
                  Diary +
                  Writing +  
                  Painting +  
                  Montage + 
                  Circus +
                  Pottery + 
                  Theater + 
                  Drawing + 
                  Dancing +  
                  Photography +
                  Genealogy + 
                  Science +
                  None +
                  No_Amateur +
                  Video_games +
                  TV +
                  Radio+
                  Library+
                  Museums +  
                  Internet + 
                  Concert, data=my_data_frame)



n <- nrow(my_data_frame)

# Indices pour la division (2/3 pour l'entraînement, 1/3 pour le test)
train_index <- sample(1:n, size = 2 * n / 3)  # 2/3 des indices pour l'entraînement

# Créer l'ensemble d'entraînement et l'ensemble de test
train_data <- my_data_frame[train_index, ]  # Enregistrement d'entraînement
test_data <- my_data_frame[-train_index, ]  # Enregistrement de test

# Vérifier les tailles
cat("Nombre d'observations dans l'ensemble d'entraînement :", nrow(train_data), "\n")
Nombre d'observations dans l'ensemble d'entraînement : 6156 
code R
cat("Nombre d'observations dans l'ensemble de test :", nrow(test_data), "\n")
Nombre d'observations dans l'ensemble de test : 3078 
code R
# Charger les bibliothèques nécessaires
library(glmnet)
library(pROC)

# Définir la matrice X pour l'entraînement
x <- model.matrix(SEXE ~ Knitting + Cards_games + Gambling + Cooking + DIY +
                  Vegetable_garden + Ornamental_garden + Fishing_hunting + 
                  Collection + Vehicle_custom + Making_music + Diary + Writing + 
                  Painting + Montage + Circus + Pottery + Theater + Drawing + 
                  Dancing + Photography + Genealogy + Science + None + No_Amateur +
                  Video_games + TV + Radio + Library + Museums + Internet + Concert,
                  train_data)[, -1]

# Variable cible y pour l'entraînement
y <- train_data$SEXE

# Ajuster le modèle LASSO avec validation croisée
cv_lasso <- cv.glmnet(x, y, alpha = 1, family = "binomial")
best_lambda <- cv_lasso$lambda.min  # Lambda optimal

# Vérifie les coefficients pour le lambda optimal
print(coef(cv_lasso, s = "lambda.min"))
33 x 1 sparse Matrix of class "dgCMatrix"
                             s1
(Intercept)       -5.521795e-01
Knitting           2.882558e+00
Cards_games        2.349254e-01
Gambling          -3.337312e-01
Cooking            1.162588e+00
DIY               -1.086545e+00
Vegetable_garden  -3.464071e-01
Ornamental_garden  3.256615e-01
Fishing_hunting   -1.450190e+00
Collection        -5.394222e-01
Vehicle_custom    -1.510231e+00
Making_music      -1.484874e-01
Diary              1.379585e+00
Writing            .           
Painting           4.704994e-01
Montage           -9.872001e-01
Circus             .           
Pottery            4.123501e-01
Theater           -8.086209e-02
Drawing           -2.584548e-02
Dancing            1.596506e+00
Photography       -4.617016e-01
Genealogy         -2.927618e-01
Science           -6.851022e-01
None              -9.890995e-02
No_Amateur         3.514798e-01
Video_games       -1.754886e-01
TV                 4.929155e-01
Radio             -1.240134e-01
Library           -9.026366e-16
Museums           -1.205780e-01
Internet           1.147213e-01
Concert            8.777883e-02
code R
print(best_lambda)
[1] 0.001109779
code R
# Préparer X_test (même traitement que pour l'entraînement)
x_test <- model.matrix(SEXE ~ Knitting + Cards_games + Gambling + Cooking + DIY +
                       Vegetable_garden + Ornamental_garden + Fishing_hunting + 
                       Collection + Vehicle_custom + Making_music + Diary + Writing + 
                       Painting + Montage + Circus + Pottery + Theater + Drawing + 
                       Dancing + Photography + Genealogy + Science + None + No_Amateur +
                       Video_games + TV + Radio + Library + Museums + Internet + Concert,
                       test_data)[, -1]

# Prédire les probabilités sur les données de test avec le meilleur lambda
preds <- predict(cv_lasso, newx = x_test, s = "lambda.min", type = "response")

# Variable cible y_test pour l'ensemble de test
y_test <- test_data$SEXE

# Calcul de l'AUC avec pROC
roc_curve <- roc(y_test, preds)
auc_value <- auc(roc_curve)
print(paste("AUC:", auc_value))
[1] "AUC: 0.873070159752027"
code R
# Calculer la courbe ROC
roc_curve <- roc(y_test, preds)

# Afficher la courbe ROC
plot(roc_curve, main = "Courbe ROC", col = "blue", lwd = 2)

# Optionnel : Ajouter la ligne de base (diagonale)
abline(a = 0, b = 1, col = "red", lty = 2)

code R
# Afficher les coefficients pour le meilleur lambda
coefficients <- coef(cv_lasso, s = "lambda.min")

# Convertir en data.frame pour une meilleure lisibilité
coefficients_df <- as.data.frame(as.matrix(coefficients))
colnames(coefficients_df) <- "Coefficient"
coefficients_df$Variable <- rownames(coefficients_df)

# Filtrer pour ne garder que les variables avec des coefficients non nuls
coefficients_non_zero <- coefficients_df[coefficients_df$Coefficient != 0, ]

# Afficher les variables gardées et leurs coefficients
print(coefficients_non_zero)
                    Coefficient          Variable
(Intercept)       -5.521795e-01       (Intercept)
Knitting           2.882558e+00          Knitting
Cards_games        2.349254e-01       Cards_games
Gambling          -3.337312e-01          Gambling
Cooking            1.162588e+00           Cooking
DIY               -1.086545e+00               DIY
Vegetable_garden  -3.464071e-01  Vegetable_garden
Ornamental_garden  3.256615e-01 Ornamental_garden
Fishing_hunting   -1.450190e+00   Fishing_hunting
Collection        -5.394222e-01        Collection
Vehicle_custom    -1.510231e+00    Vehicle_custom
Making_music      -1.484874e-01      Making_music
Diary              1.379585e+00             Diary
Painting           4.704994e-01          Painting
Montage           -9.872001e-01           Montage
Pottery            4.123501e-01           Pottery
Theater           -8.086209e-02           Theater
Drawing           -2.584548e-02           Drawing
Dancing            1.596506e+00           Dancing
Photography       -4.617016e-01       Photography
Genealogy         -2.927618e-01         Genealogy
Science           -6.851022e-01           Science
None              -9.890995e-02              None
No_Amateur         3.514798e-01        No_Amateur
Video_games       -1.754886e-01       Video_games
TV                 4.929155e-01                TV
Radio             -1.240134e-01             Radio
Library           -9.026366e-16           Library
Museums           -1.205780e-01           Museums
Internet           1.147213e-01          Internet
Concert            8.777883e-02           Concert
code R
###CREATION INDICE

# Récupérer les coefficients du modèle pour le meilleur lambda
coefficients <- coef(cv_lasso, s = "lambda.min")

# Convertir en data.frame et filtrer les variables non nulles
coefficients_df <- as.data.frame(as.matrix(coefficients))
colnames(coefficients_df) <- "Coefficient"
coefficients_df$Variable <- rownames(coefficients_df)
coefficients_non_zero <- coefficients_df[coefficients_df$Coefficient != 0, ]

# Préparer la matrice X des pratiques pour l'ensemble de test (ou d'entraînement si nécessaire)
x_test <- model.matrix(SEXE ~ Knitting + Cards_games + Gambling + Cooking + DIY +
                       Vegetable_garden + Ornamental_garden + Fishing_hunting + 
                       Collection + Vehicle_custom + Making_music + Diary + Writing + 
                       Painting + Montage + Circus + Pottery + Theater + Drawing + 
                       Dancing + Photography + Genealogy + Science + None + No_Amateur +
                       Video_games + TV + Radio + Library + Museums + Internet + Concert,
                       test_data)[, -1]  # Assurez-vous d'enlever l'intercept avec [, -1]

# Sélectionner les variables pertinentes
matched_vars <- intersect(rownames(coefficients_non_zero), colnames(x_test))
x_test <- x_test[, matched_vars, drop = FALSE]
coef_vector <- coefficients_non_zero$Coefficient[matched_vars]  # Associer les coefficients

print(paste("Dimension de x_test :", dim(x_test)[1], "x", dim(x_test)[2]))
[1] "Dimension de x_test : 3078 x 30"
code R
print(paste("Longueur de coef_vector :", length(coef_vector)))
[1] "Longueur de coef_vector : 30"
code R
print(paste("Nombre de valeurs NA dans x_test :", sum(is.na(x_test))))
[1] "Nombre de valeurs NA dans x_test : 0"
code R
print(paste("Nombre de valeurs NA dans coef_vector :", sum(is.na(coef_vector))))
[1] "Nombre de valeurs NA dans coef_vector : 30"
code R
print("Variables dans coefficients_non_zero :")
[1] "Variables dans coefficients_non_zero :"
code R
print(coefficients_non_zero$Variable)
 [1] "(Intercept)"       "Knitting"          "Cards_games"      
 [4] "Gambling"          "Cooking"           "DIY"              
 [7] "Vegetable_garden"  "Ornamental_garden" "Fishing_hunting"  
[10] "Collection"        "Vehicle_custom"    "Making_music"     
[13] "Diary"             "Painting"          "Montage"          
[16] "Pottery"           "Theater"           "Drawing"          
[19] "Dancing"           "Photography"       "Genealogy"        
[22] "Science"           "None"              "No_Amateur"       
[25] "Video_games"       "TV"                "Radio"            
[28] "Library"           "Museums"           "Internet"         
[31] "Concert"          
code R
print("Colonnes de x_test :")
[1] "Colonnes de x_test :"
code R
print(colnames(x_test))
 [1] "Knitting"          "Cards_games"       "Gambling"         
 [4] "Cooking"           "DIY"               "Vegetable_garden" 
 [7] "Ornamental_garden" "Fishing_hunting"   "Collection"       
[10] "Vehicle_custom"    "Making_music"      "Diary"            
[13] "Painting"          "Montage"           "Pottery"          
[16] "Theater"           "Drawing"           "Dancing"          
[19] "Photography"       "Genealogy"         "Science"          
[22] "None"              "No_Amateur"        "Video_games"      
[25] "TV"                "Radio"             "Library"          
[28] "Museums"           "Internet"          "Concert"          
code R
# Exclure "(Intercept)" des coefficients
coefficients_filtered <- coefficients_non_zero[coefficients_non_zero$Variable != "(Intercept)", ]

# Aligner les coefficients sur x_test
coef_vector <- coefficients_filtered$Coefficient[match(colnames(x_test), coefficients_filtered$Variable)]



# Calcul des scores
test_data$score <- x_test %*% coef_vector

# Normalisation des scores
min_score <- min(test_data$score, na.rm = TRUE)
max_score <- max(test_data$score, na.rm = TRUE)

if (max_score > min_score) {
  test_data$score_normalise <- (test_data$score - min_score) / (max_score - min_score)
} else {
  test_data$score_normalise <- 0
}

# Visualisation avec ggplot2
library(ggplot2)
ggplot(test_data, aes(x = score_normalise, color = SEXE, fill = SEXE)) +
  geom_density(alpha = 0.4) + 
  labs(title = "Densité du Score Normalisé", x = "Score Normalisé", y = "Densité") +
  scale_fill_manual(values = c("blue", "pink")) +
  scale_color_manual(values = c("blue", "pink")) +
  theme_minimal()

code R
# Préparer la matrice X pour l'ensemble complet (train + test)
x_full <- model.matrix(SEXE ~ Knitting + Cards_games + Gambling + Cooking + DIY +
                       Vegetable_garden + Ornamental_garden + Fishing_hunting + 
                       Collection + Vehicle_custom + Making_music + Diary + Writing + 
                       Painting + Montage + Circus + Pottery + Theater + Drawing + 
                       Dancing + Photography + Genealogy + Science + None + No_Amateur +
                       Video_games + TV + Radio + Library + Museums + Internet + Concert,
                       my_data_frame)[, -1]  # Exclure l'intercept

# Sélectionner les variables pertinentes (identiques à celles du modèle LASSO)
matched_vars <- intersect(rownames(coefficients_non_zero), colnames(x_full))
x_full <- x_full[, matched_vars, drop = FALSE]  # Garder uniquement les variables sélectionnées
# Aligner les coefficients sur x_full
coef_vector_full <- coefficients_filtered$Coefficient[match(colnames(x_full), coefficients_filtered$Variable)]

# Calcul des scores pour l'ensemble complet
my_data_frame$score <- x_full %*% coef_vector_full
# Normalisation des scores
min_score_full <- min(my_data_frame$score, na.rm = TRUE)
max_score_full <- max(my_data_frame$score, na.rm = TRUE)

if (max_score_full > min_score_full) {
  my_data_frame$score_normalise <- (my_data_frame$score - min_score_full) / (max_score_full - min_score_full)
} else {
  my_data_frame$score_normalise <- 0
}
# Visualisation du score normalisé pour l'ensemble complet avec ggplot2
library(ggplot2)
ggplot(my_data_frame, aes(x = score_normalise, color = SEXE, fill = SEXE)) +
  geom_density(alpha = 0.4) + 
  labs(title = "Densité du Score Normalisé (Ensemble Complet)", x = "Score Normalisé", y = "Densité") +
  scale_fill_manual(values = c("blue", "pink")) +
  scale_color_manual(values = c("blue", "pink")) +
  theme_minimal()

1.2 Comparaisons

code R
# Calculer la corrélation entre indice_normalise et score_normalise
correlation_value <- cor(my_data_frame$identity, my_data_frame$score_normalise, use = "complete.obs")
print(paste("Corrélation entre indice_normalise et score_normalise :", correlation_value))
[1] "Corrélation entre indice_normalise et score_normalise : -0.645386557366119"
code R
# Visualiser la relation entre indice_normalise et score_normalise
library(ggplot2)
ggplot(my_data_frame, aes(x = identity, y = score_normalise)) +
  geom_point(alpha = 0.5) +
  labs(title = "Relation entre indice_normalise et score_normalise", 
       x = "Indice Normalisé", 
       y = "Score Normalisé") +
  theme_minimal() +
  geom_smooth(method = "lm", col = "red", se = FALSE)  # Ajouter une droite de régression linéaire

code R
n_2 <- nrow(my_data_frame)

# Indices pour la division (2/3 pour l'entraînement, 1/3 pour le test)
train_index <- sample(1:n_2, size = 2 * n_2 / 3)  # 2/3 des indices pour l'entraînement

# Créer l'ensemble d'entraînement et l'ensemble de test
train_data <- my_data_frame[train_index, ]  # Enregistrement d'entraînement
test_data <- my_data_frame[-train_index, ]  # Enregistrement de test

# Vérifier les tailles
cat("Nombre d'observations dans l'ensemble d'entraînement :", nrow(train_data), "\n")
Nombre d'observations dans l'ensemble d'entraînement : 6156 
code R
cat("Nombre d'observations dans l'ensemble de test :", nrow(test_data), "\n")
Nombre d'observations dans l'ensemble de test : 3078 
code R
# Charger les bibliothèques nécessaires
library(glmnet)
library(pROC)

# Définir la matrice X pour l'entraînement
x_2 <- model.matrix(SEXE ~ Knitting + Cards_games + Gambling + Cooking + DIY +
                  Vegetable_garden + Fishing_hunting + 
                  Collection + Vehicle_custom + Making_music + Diary + Writing + 
                  Painting + Montage +  Pottery + Theater + Drawing + 
                  Dancing + Photography + Genealogy + Science + None + 
                  Video_games + Library +  Concert,
                  train_data)[, -1]

# Variable cible y pour l'entraînement
y_2 <- train_data$SEXE

# Ajuster le modèle LASSO avec validation croisée
cv_lasso_2 <- cv.glmnet(x_2, y_2, alpha = 1, family = "binomial")
best_lambda_2 <- cv_lasso_2$lambda.min  # Lambda optimal

# Vérifie les coefficients pour le lambda optimal
print(coef(cv_lasso_2, s = "lambda.min"))
26 x 1 sparse Matrix of class "dgCMatrix"
                          s1
(Intercept)       0.11260864
Knitting          3.01800722
Cards_games       0.16677278
Gambling         -0.34566529
Cooking           1.12798392
DIY              -1.01219567
Vegetable_garden -0.30731087
Fishing_hunting  -1.40787072
Collection       -0.69110001
Vehicle_custom   -1.68777530
Making_music     -0.13299836
Diary             1.46711705
Writing          -0.16958639
Painting          0.48089960
Montage          -1.00019427
Pottery           0.52939803
Theater          -0.06515215
Drawing          -0.15971724
Dancing           1.52751004
Photography      -0.50794702
Genealogy        -0.33199701
Science          -0.87343121
None             -0.11072248
Video_games      -0.16572230
Library          -0.01094454
Concert           0.01127190
code R
print(best_lambda_2)
[1] 0.0005315669
code R
# Préparer X_test (même traitement que pour l'entraînement)
x_test_2 <- model.matrix(SEXE ~ Knitting + Cards_games + Gambling + Cooking + DIY +
                  Vegetable_garden + Fishing_hunting + 
                  Collection + Vehicle_custom + Making_music + Diary + Writing + 
                  Painting + Montage +  Pottery + Theater + Drawing + 
                  Dancing + Photography + Genealogy + Science + None + 
                  Video_games + Library +  Concert,
                       test_data)[, -1]

# Prédire les probabilités sur les données de test avec le meilleur lambda
preds_2 <- predict(cv_lasso_2, newx = x_test_2, s = "lambda.min", type = "response")

# Variable cible y_test pour l'ensemble de test
y_test_2 <- test_data$SEXE

# Calcul de l'AUC avec pROC
roc_curve_2 <- roc(y_test_2, preds_2)
auc_value_2 <- auc(roc_curve_2)
print(paste("AUC:", auc_value_2))
[1] "AUC: 0.868197882559135"
code R
# Calculer la courbe ROC
roc_curve_2 <- roc(y_test_2, preds_2)

# Afficher la courbe ROC
plot(roc_curve_2, main = "Courbe ROC 2", col = "blue", lwd = 2)

# Optionnel : Ajouter la ligne de base (diagonale)
abline(a = 0, b = 1, col = "red", lty = 2)

code R
# Afficher les coefficients pour le meilleur lambda
coefficients_2 <- coef(cv_lasso_2, s = "lambda.min")

# Convertir en data.frame pour une meilleure lisibilité
coefficients_df_2 <- as.data.frame(as.matrix(coefficients_2))
colnames(coefficients_df_2) <- "Coefficient"
coefficients_df_2$Variable <- rownames(coefficients_df_2)

# Filtrer pour ne garder que les variables avec des coefficients non nuls
coefficients_non_zero_2 <- coefficients_df_2[coefficients_df_2$Coefficient != 0, ]

# Afficher les variables gardées et leurs coefficients
print(coefficients_non_zero_2)
                 Coefficient         Variable
(Intercept)       0.11260864      (Intercept)
Knitting          3.01800722         Knitting
Cards_games       0.16677278      Cards_games
Gambling         -0.34566529         Gambling
Cooking           1.12798392          Cooking
DIY              -1.01219567              DIY
Vegetable_garden -0.30731087 Vegetable_garden
Fishing_hunting  -1.40787072  Fishing_hunting
Collection       -0.69110001       Collection
Vehicle_custom   -1.68777530   Vehicle_custom
Making_music     -0.13299836     Making_music
Diary             1.46711705            Diary
Writing          -0.16958639          Writing
Painting          0.48089960         Painting
Montage          -1.00019427          Montage
Pottery           0.52939803          Pottery
Theater          -0.06515215          Theater
Drawing          -0.15971724          Drawing
Dancing           1.52751004          Dancing
Photography      -0.50794702      Photography
Genealogy        -0.33199701        Genealogy
Science          -0.87343121          Science
None             -0.11072248             None
Video_games      -0.16572230      Video_games
Library          -0.01094454          Library
Concert           0.01127190          Concert
code R
###CREATION INDICE

# Récupérer les coefficients du modèle pour le meilleur lambda
coefficients_2 <- coef(cv_lasso_2, s = "lambda.min")

# Convertir en data.frame et filtrer les variables non nulles
coefficients_df_2 <- as.data.frame(as.matrix(coefficients_2))
colnames(coefficients_df_2) <- "Coefficient"
coefficients_df_2$Variable <- rownames(coefficients_df_2)
coefficients_non_zero_2 <- coefficients_df_2[coefficients_df_2$Coefficient != 0, ]

# Préparer la matrice X des pratiques pour l'ensemble de test (ou d'entraînement si nécessaire)
x_test_2 <- model.matrix(SEXE ~ Knitting + Cards_games + Gambling + Cooking + DIY +
                  Vegetable_garden + Fishing_hunting + 
                  Collection + Vehicle_custom + Making_music + Diary + Writing + 
                  Painting + Montage +  Pottery + Theater + Drawing + 
                  Dancing + Photography + Genealogy + Science + None + 
                  Video_games + Library +  Concert,
                       test_data)[, -1]  # Assurez-vous d'enlever l'intercept avec [, -1]

# Sélectionner les variables pertinentes
matched_vars_2 <- intersect(rownames(coefficients_non_zero_2), colnames(x_test_2))
x_test_2 <- x_test_2[, matched_vars_2, drop = FALSE]
coef_vector_2 <- coefficients_non_zero_2$Coefficient[matched_vars_2]  # Associer les coefficients

#print(paste("Dimension de x_test :", dim(x_test)[1], "x", dim(x_test)[2]))
#print(paste("Longueur de coef_vector :", length(coef_vector)))
#print(paste("Nombre de valeurs NA dans x_test :", sum(is.na(x_test))))
#print(paste("Nombre de valeurs NA dans coef_vector :", sum(is.na(coef_vector))))

#print("Variables dans coefficients_non_zero :")
#print(coefficients_non_zero$Variable)

#print("Colonnes de x_test :")
#print(colnames(x_test))


# Exclure "(Intercept)" des coefficients
coefficients_filtered_2 <- coefficients_non_zero_2[coefficients_non_zero_2$Variable != "(Intercept)", ]

# Aligner les coefficients sur x_test
coef_vector_2 <- coefficients_filtered_2$Coefficient[match(colnames(x_test_2), coefficients_filtered_2$Variable)]



# Calcul des scores
test_data$score_2 <- x_test_2 %*% coef_vector_2

# Normalisation des scores
min_score_2 <- min(test_data$score_2, na.rm = TRUE)
max_score_2 <- max(test_data$score_2, na.rm = TRUE)

if (max_score_2 > min_score_2) {
  test_data$score_normalise_2 <- (test_data$score_2 - min_score_2) / (max_score_2 - min_score_2)
} else {
  test_data$score_normalise_2 <- 0
}

# Visualisation avec ggplot2
library(ggplot2)
ggplot(test_data, aes(x = score_normalise_2, color = SEXE, fill = SEXE)) +
  geom_density(alpha = 0.4) + 
  labs(title = "Densité du Score Normalisé 2", x = "Score Normalisé 2", y = "Densité") +
  scale_fill_manual(values = c("blue", "pink")) +
  scale_color_manual(values = c("blue", "pink")) +
  theme_minimal()

code R
# Préparer la matrice X pour l'ensemble complet (train + test)
x_full_2 <- model.matrix(SEXE ~ Knitting + Cards_games + Gambling + Cooking + DIY +
                  Vegetable_garden + Fishing_hunting + 
                  Collection + Vehicle_custom + Making_music + Diary + Writing + 
                  Painting + Montage +  Pottery + Theater + Drawing + 
                  Dancing + Photography + Genealogy + Science + None + 
                  Video_games + Library +  Concert,
                       my_data_frame)[, -1]  # Exclure l'intercept

# Sélectionner les variables pertinentes (identiques à celles du modèle LASSO)
matched_vars_2 <- intersect(rownames(coefficients_non_zero_2), colnames(x_full_2))
x_full_2 <- x_full_2[, matched_vars_2, drop = FALSE]  # Garder uniquement les variables sélectionnées
# Aligner les coefficients sur x_full
coef_vector_full_2 <- coefficients_filtered_2$Coefficient[match(colnames(x_full_2), coefficients_filtered_2$Variable)]

# Calcul des scores pour l'ensemble complet
my_data_frame$score_2 <- x_full_2 %*% coef_vector_full_2
# Normalisation des scores
min_score_full_2 <- min(my_data_frame$score_2, na.rm = TRUE)
max_score_full_2 <- max(my_data_frame$score_2, na.rm = TRUE)

if (max_score_full_2 > min_score_full_2) {
  my_data_frame$score_normalise_2 <- (my_data_frame$score_2 - min_score_full_2) / (max_score_full_2 - min_score_full_2)
} else {
  my_data_frame$score_normalise_2 <- 0
}
# Visualisation du score normalisé pour l'ensemble complet avec ggplot2
library(ggplot2)
ggplot(my_data_frame, aes(x = score_normalise_2, color = SEXE, fill = SEXE)) +
  geom_density(alpha = 0.4) + 
  labs(title = "Densité du Score Normalisé 2(Ensemble Complet)", x = "Score Normalisé", y = "Densité") +
  scale_fill_manual(values = c("blue", "pink")) +
  scale_color_manual(values = c("blue", "pink")) +
  theme_minimal()

code R
# Calculer la corrélation entre indice_normalise et score_normalise
correlation_value_2 <- cor(my_data_frame$identity, my_data_frame$score_normalise_2, use = "complete.obs")
print(paste("Corrélation entre indice_normalise et score_normalise 2:", correlation_value))
[1] "Corrélation entre indice_normalise et score_normalise 2: -0.645386557366119"
code R
# Visualiser la relation entre indice_normalise et score_normalise
library(ggplot2)
ggplot(my_data_frame, aes(x = identity, y = score_normalise_2)) +
  geom_point(alpha = 0.5) +
  labs(title = "Relation entre indice_normalise et score_normalise 2", 
       x = "Indice Normalisé", 
       y = "Score Normalisé") +
  theme_minimal() +
  geom_smooth(method = "lm", col = "red", se = FALSE)  # Ajouter une droite de régression linéaire

code R
# Créer un data frame pour stocker les résultats
result_table_score_2 <- data.frame(Activity_2 = character(0), Accuracy_2 = numeric(0))

cultural_activities_2 <- c(
   "Knitting", "Cards_games", "Gambling", "Cooking", "DIY",
   "Vegetable_garden", "Fishing_hunting", "Collection", "Vehicle_custom",
   "Making_music", "Diary", "Writing", "Painting", "Montage",
   "Pottery", "Theater", "Drawing", "Dancing", "Photography",
   "Genealogy", "Science", "None", "Video_games", "Library", "Concert"
)

for (activity_2 in cultural_activities_2) {
  # Effectuer une régression Probit pour l'activité culturelle actuelle
  model_formula_2 <- as.formula(paste(activity_2, "~ score_normalise_2"))
  model_2 <- glm(model_formula_2, data = my_data_frame, family = binomial(link = "probit"))
 
  # Calculer les prédictions
  predicted_2 <- ifelse(predict(model_2, type = "response") >= 0.5, 1, 0)
 
  # Calculer la précision
  correct_predictions_2 <- sum(predicted_2 == my_data_frame[[activity_2]])
  total_predictions_2 <- nrow(my_data_frame)
  accuracy_2 <- (correct_predictions_2 / total_predictions_2) * 100
 
  # Ajouter le résultat au data frame
  result_table_score_2 <- rbind(result_table_score_2, data.frame(Activity_2 = activity_2, Accuracy_2 = accuracy_2))
}

# Afficher le tableau des résultats
print(result_table_score_2)
         Activity_2 Accuracy_2
1          Knitting   92.92831
2       Cards_games   56.44358
3          Gambling   78.77410
4           Cooking   68.72428
5               DIY   59.97401
6  Vegetable_garden   72.07061
7   Fishing_hunting   90.25341
8        Collection   92.67923
9    Vehicle_custom   96.58869
10     Making_music   65.97358
11            Diary   84.53541
12          Writing   87.48105
13         Painting   78.68746
14          Montage   84.54624
15          Pottery   89.50617
16          Theater   86.10570
17          Drawing   76.72731
18          Dancing   76.80312
19      Photography   74.53974
20        Genealogy   88.21746
21          Science   88.30409
22             None   71.20425
23      Video_games   61.47932
24          Library   83.47412
25          Concert   82.06628
code R
my_data_frame$score_scale <- cut(
  my_data_frame$score_normalise,
  breaks = 7,  # Automatically divide into 7 slices
  labels = c("Very Masculine", "1", "2", "3", "4", "5", "Very Feminine"), 
  include.lowest = TRUE  
)
table1 <-
  my_data_frame |> 
  tbl_summary(include = c(score_scale),
  by=Sex ,)|> 
    add_p()
table1
Characteristic Men
N = 4,1621
Women
N = 5,0721
p-value2
score_scale

<0.001
    Very Masculine 18 (0.4%) 0 (0%)
    1 405 (9.7%) 15 (0.3%)
    2 1,871 (45%) 412 (8.1%)
    3 1,659 (40%) 2,075 (41%)
    4 189 (4.5%) 1,646 (32%)
    5 19 (0.5%) 792 (16%)
    Very Feminine 1 (<0.1%) 132 (2.6%)
1 n (%)
2 Pearson’s Chi-squared test
code R
# Proportions de 'score_scale' par genre
table_score_gender <- table(my_data_frame$score_scale, my_data_frame$Sex)
table_score_gender_percent <- prop.table(table_score_gender, 2) * 100  # Calcul par genre
table_score_gender_percent
                
                         Men       Women
  Very Masculine  0.43248438  0.00000000
  1               9.73089861  0.29574132
  2              44.95434887  8.12302839
  3              39.86064392 40.91088328
  4               4.54108602 32.45268139
  5               0.45651129 15.61514196
  Very Feminine   0.02402691  2.60252366
code R
# Proportions de 'satisfaction' par genre
table_satisfaction_gender <- table(my_data_frame$satisfaction, my_data_frame$Sex)
table_satisfaction_gender_percent <- prop.table(table_satisfaction_gender, 2) * 100  # Calcul par genre
table_satisfaction_gender_percent
        
              Men    Women
  High   35.15399 32.66917
  Low    35.80366 39.75143
  Medium 29.04235 27.57940
code R
# Charger les bibliothèques nécessaires
library(ggplot2)
library(plotly)
library(scales)

# Création du graphique avec ggplot
p <- ggplot(as.data.frame(table_score_gender_percent), 
            aes(x = Var1, y = Freq, fill = Var2, text = paste0("Proportion: ", percent(Freq / 100)))) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(values = c("blue", "pink")) +
  labs(title = "Répartition des scores par genre",
       x = "Score (1-7)", y = "Proportion (%)") +
  scale_y_continuous(labels = percent_format(scale = 1)) +
  theme_minimal() +
  theme(legend.title = element_blank())

# Convertir en graphique interactif avec ggplotly
ggplotly(p, tooltip = "text")
code R
table1 <-
  my_data_frame |> 
  tbl_summary(include = c(identity_scale),
  by=Sex ,)|> 
    add_p()
table1
Characteristic Men
N = 4,1621
Women
N = 5,0721
p-value2
identity_scale

<0.001
    Very Feminine 9 (0.2%) 218 (4.3%)
    1 403 (9.7%) 1,506 (30%)
    2 2,139 (51%) 2,669 (53%)
    3 985 (24%) 576 (11%)
    4 509 (12%) 90 (1.8%)
    5 96 (2.3%) 10 (0.2%)
    Very Masculine 21 (0.5%) 3 (<0.1%)
1 n (%)
2 Pearson’s Chi-squared test
code R
# Proportions de 'score_scale' par genre
table_identity_gender <- table(my_data_frame$identity_scale, my_data_frame$Sex)
table_identity_gender_percent <- prop.table(table_identity_gender, 2) * 100  # Calcul par genre
table_identity_gender_percent
                
                         Men       Women
  Very Feminine   0.21624219  4.29810726
  1               9.68284479 29.69242902
  2              51.39356079 52.62223975
  3              23.66650649 11.35646688
  4              12.22969726  1.77444795
  5               2.30658337  0.19716088
  Very Masculine  0.50456511  0.05914826
code R
# Création du graphique avec ggplot
p <- ggplot(as.data.frame(table_identity_gender_percent), 
            aes(x = Var1, y = Freq, fill = Var2, text = paste0("Proportion: ", percent(Freq / 100)))) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(values = c("blue", "pink")) +
  labs(title = "Répartition de l'identity scale par genre",
       x = "Score (1-7)", y = "Proportion (%)") +
  scale_y_continuous(labels = percent_format(scale = 1)) +
  theme_minimal() +
  theme(legend.title = element_blank())

# Convertir en graphique interactif avec ggplotly
ggplotly(p, tooltip = "text")
code R
table1 <-
  my_data_frame |> 
  tbl_summary(include = c(CLASSIF),
  by=score_scale ,)|> 
    add_p()
table1
Characteristic Very Masculine
N = 181
1
N = 4201
2
N = 2,2831
3
N = 3,7341
4
N = 1,8351
5
N = 8111
Very Feminine
N = 1331
p-value
CLASSIF







    1 0 (0%) 14 (8.2%) 72 (9.1%) 70 (6.4%) 23 (4.3%) 6 (3.2%) 1 (3.1%)
    2 3 (38%) 63 (37%) 236 (30%) 221 (20%) 59 (11%) 14 (7.5%) 1 (3.1%)
    3 0 (0%) 24 (14%) 105 (13%) 99 (9.1%) 39 (7.4%) 13 (7.0%) 2 (6.3%)
    4 3 (38%) 22 (13%) 83 (10%) 97 (8.9%) 56 (11%) 17 (9.1%) 5 (16%)
    5 0 (0%) 26 (15%) 140 (18%) 191 (18%) 103 (19%) 44 (24%) 10 (31%)
    6 1 (13%) 20 (12%) 151 (19%) 380 (35%) 245 (46%) 87 (47%) 13 (41%)
    7 0 (0%) 0 (0%) 4 (0.5%) 12 (1.1%) 2 (0.4%) 2 (1.1%) 0 (0%)
    8 1 (13%) 1 (0.6%) 4 (0.5%) 16 (1.5%) 3 (0.6%) 4 (2.1%) 0 (0%)
    9 0 (0%) 0 (0%) 0 (0%) 2 (0.2%) 0 (0%) 0 (0%) 0 (0%)
    Unknown 10 250 1,488 2,646 1,305 624 101
1 n (%)
code R
# Proportions de 'score_scale' par profession
table_identity_gender <- table(my_data_frame$CLASSIF, my_data_frame$score_scale)
table_identity_gender_percent <- prop.table(table_identity_gender, 2) * 100  # Calcul par genre
table_identity_gender_percent
   
    Very Masculine          1          2          3          4          5
  1      0.0000000  8.2352941  9.0566038  6.4338235  4.3396226  3.2085561
  2     37.5000000 37.0588235 29.6855346 20.3125000 11.1320755  7.4866310
  3      0.0000000 14.1176471 13.2075472  9.0992647  7.3584906  6.9518717
  4     37.5000000 12.9411765 10.4402516  8.9154412 10.5660377  9.0909091
  5      0.0000000 15.2941176 17.6100629 17.5551471 19.4339623 23.5294118
  6     12.5000000 11.7647059 18.9937107 34.9264706 46.2264151 46.5240642
  7      0.0000000  0.0000000  0.5031447  1.1029412  0.3773585  1.0695187
  8     12.5000000  0.5882353  0.5031447  1.4705882  0.5660377  2.1390374
  9      0.0000000  0.0000000  0.0000000  0.1838235  0.0000000  0.0000000
   
    Very Feminine
  1     3.1250000
  2     3.1250000
  3     6.2500000
  4    15.6250000
  5    31.2500000
  6    40.6250000
  7     0.0000000
  8     0.0000000
  9     0.0000000
code R
# Création du graphique avec ggplot
p <- ggplot(as.data.frame(table_identity_gender_percent), 
            aes(x = Var1, y = Freq, fill = Var2, text = paste0("Proportion: ", percent(Freq / 100)))) +
  geom_bar(stat = "identity", position = "dodge") +
 
  labs(title = "Répartition des professions par score scale",
       x = "Profession", y = "Proportion (%)") +
  scale_y_continuous(labels = percent_format(scale = 1)) +
  theme_minimal() +
  theme(legend.title = element_blank())

# Convertir en graphique interactif avec ggplotly
ggplotly(p, tooltip = "text")
code R
library(ggplot2)
library(dplyr)
library(plotly)

# Créer un tableau de fréquence pour la heatmap
heatmap_data <- my_data_frame %>%
  count(score_scale, satisfaction)  # Compte les occurrences

# Ajouter une colonne pour l'affichage interactif
heatmap_data <- heatmap_data %>%
  mutate(info_text = paste0("Score: ", score_scale, 
                            "\nSatisfaction: ", satisfaction, 
                            "\nNombre d'observations: ", n))

# Créer la heatmap avec ggplot
p <- ggplot(heatmap_data, aes(x = as.factor(score_scale), 
                              y = as.factor(satisfaction), 
                              fill = n, text = info_text)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "blue") +
  labs(title = "Heatmap de la relation entre satisfaction et score",
       x = "Score Scale", y = "Satisfaction", fill = "Fréquence") +
  theme_minimal()

# Rendre le graphique interactif avec des infobulles personnalisées
ggplotly(p, tooltip = "text")
code R
# Créer un tableau de fréquence pour la heatmap
heatmap_data <- my_data_frame %>%
  count(identity_scale, satisfaction)  # Compte les occurrences

# Ajouter une colonne pour l'affichage interactif
heatmap_data <- heatmap_data %>%
  mutate(info_text = paste0("Identity: ", identity_scale, 
                            "\nSatisfaction: ", satisfaction, 
                            "\nNombre d'observations: ", n))

# Créer la heatmap avec ggplot
p <- ggplot(heatmap_data, aes(x = as.factor(identity_scale), 
                              y = as.factor(satisfaction), 
                              fill = n, text = info_text)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "blue") +
  labs(title = "Heatmap de la relation entre satisfaction et identity scale",
       x = "Score Scale", y = "Satisfaction", fill = "Fréquence") +
  theme_minimal()

# Rendre le graphique interactif avec des infobulles personnalisées
ggplotly(p, tooltip = "text")

1.3 Ecarts aux normes et satisfaction

code R
mean_gender <- tapply(my_data_frame$identity, my_data_frame$Sex, mean)
sd_gender <- tapply(my_data_frame$identity, my_data_frame$Sex, sd)

my_data_frame$distance_abs<- abs((my_data_frame$identity - mean_gender[my_data_frame$Sex]) / sd_gender[my_data_frame$Sex])

ggplot(my_data_frame, aes(x = distance_abs, color = Sex, fill = Sex)) +
  geom_density(alpha = 0.4) +
  labs(title = "Density of Distance to the Norm, by Sex", 
       x = "Distance to the Norm (Z-score)", y = "Density") +
  scale_fill_manual(values = c("blue", "pink")) +
  scale_color_manual(values = c("blue", "pink")) +
  theme_minimal()

code R
# Charger dplyr
library(dplyr)

# Assurer que Sex est un facteur avec les niveaux appropriés
my_data_frame$Sex <- factor(my_data_frame$Sex, levels = c("Men", "Women"))

# Calculer les moyennes et écarts-types par sexe avec dplyr
mean_sd_by_sex <- my_data_frame %>%
  group_by(Sex) %>%
  summarise(
    mean_gender = mean(score_normalise, na.rm = TRUE),
    sd_gender = sd(score_normalise, na.rm = TRUE)
  )

# Joindre les moyennes et écarts-types au dataframe original
my_data_frame <- my_data_frame %>%
  left_join(mean_sd_by_sex, by = "Sex")

# Calculer la distance absolue à la moyenne (Z-score)
my_data_frame <- my_data_frame %>%
  mutate(distance_abs_2 = abs((score_normalise - mean_gender) / sd_gender))

# Visualisation
ggplot(my_data_frame, aes(x = distance_abs_2, color = Sex, fill = Sex)) +
  geom_density(alpha = 0.4) +
  labs(title = "Density of Distance to the Norm (Score), by Sex", 
       x = "Distance to the Norm (Z-score)", y = "Density") +
  scale_fill_manual(values = c("blue", "pink")) +
  scale_color_manual(values = c("blue", "pink")) +
  theme_minimal() +
  geom_vline(data = my_data_frame, aes(xintercept = mean(distance_abs_2[Sex == "Men"]), color = "Men"), linetype = "dashed") +
  geom_vline(data = my_data_frame, aes(xintercept = mean(distance_abs_2[Sex == "Women"]), color = "Women"), linetype = "dashed") +
  theme(legend.title = element_blank())  # Enlever le titre de la légende

code R
# Charger le package MASS pour polr
library(MASS)

# Modèle de régression ordinale (polr)
model <- polr(as.factor(satisfaction) ~ distance_abs_2, data = my_data_frame, method = "logistic")

# Résumé du modèle
summary(model)
Call:
polr(formula = as.factor(satisfaction) ~ distance_abs_2, data = my_data_frame, 
    method = "logistic")

Coefficients:
                 Value Std. Error t value
distance_abs_2 0.02273    0.03199  0.7105

Intercepts:
           Value    Std. Error t value 
High|Low    -0.6543   0.0340   -19.2524
Low|Medium   0.9511   0.0348    27.3139

Residual Deviance: 20135.53 
AIC: 20141.53 
(9 observations effacées parce que manquantes)
code R
# Modèle de régression ordinale avec interaction entre sexe et distance aux normes
model_interaction <- polr(as.factor(satisfaction) ~ Sex * distance_abs_2, data = my_data_frame, method = "logistic")

# Résumé du modèle
summary(model_interaction)
Call:
polr(formula = as.factor(satisfaction) ~ Sex * distance_abs_2, 
    data = my_data_frame, method = "logistic")

Coefficients:
                           Value Std. Error t value
SexWomen                 0.11297    0.06463   1.748
distance_abs_2           0.07744    0.04621   1.676
SexWomen:distance_abs_2 -0.10629    0.06407  -1.659

Intercepts:
           Value    Std. Error t value 
High|Low    -0.5961   0.0477   -12.4889
Low|Medium   1.0098   0.0486    20.7948

Residual Deviance: 20132.29 
AIC: 20142.29 
(9 observations effacées parce que manquantes)
code R
# Charger la bibliothèque Plotly
library(plotly)

# Regroupement de distance_abs_2 en catégories (bins)
my_data_frame$distance_abs_2_bins <- cut(my_data_frame$distance_abs_2, 
                                          breaks = 10,  # Divise en 10 intervalles
                                          labels = FALSE, include.lowest = TRUE)

# Comptage des occurrences pour chaque combinaison de Sexe, Satisfaction et Distance
library(dplyr)
count_data <- my_data_frame %>%
  count(Sex, satisfaction, distance_abs_2_bins)

# Création de la heatmap interactive avec Plotly
fig <- plot_ly(count_data, 
               x = ~distance_abs_2_bins, 
               y = ~as.factor(satisfaction), 
               z = ~n, 
               type = "heatmap", 
               colors = colorRamp(c("white", "blue")), 
               colorbar = list(title = "Nombre d'individus"))

# Ajouter des titres
fig <- fig %>%
  layout(title = "Heatmap de Satisfaction par Sexe et Distance aux Normes",
         xaxis = list(title = "Distance aux Normes (Binned)"),
         yaxis = list(title = "Satisfaction"))

# Afficher la heatmap interactive
fig
code R
ggplot(my_data_frame, aes(x = distance_abs_2, y = as.factor(satisfaction), color = Sex)) +
  geom_point(alpha = 0.7) +
  labs(title = "Interaction entre Sexe et Distance aux Normes sur Satisfaction",
       x = "Distance aux Normes (Z-score)", y = "Satisfaction") +
  scale_color_manual(values = c("blue", "pink")) +
  theme_minimal() +
  facet_wrap(~Sex)

code R
# 1. Charger les librairies nécessaires
library(MASS)
library(ggplot2)

# 2. Vérifier les valeurs manquantes dans le dataframe
#summary(my_data_frame)

# 3. Supprimer les lignes contenant des valeurs manquantes dans les colonnes pertinentes
my_data_frame_clean <- my_data_frame[complete.cases(my_data_frame[c("satisfaction", "Sex", "distance_abs_2")]), ]

# Vérification de la taille du dataframe nettoyé
nrow(my_data_frame_clean)  # Assure-toi que la taille est correcte
[1] 9225
code R
# 4. Convertir les variables catégoriques en facteurs avec les niveaux appropriés
my_data_frame_clean$satisfaction <- factor(my_data_frame_clean$satisfaction, levels = c("Low", "Medium", "High"))
my_data_frame_clean$Sex <- factor(my_data_frame_clean$Sex, levels = c("Men", "Women"))

# 5. Ajuster le modèle de régression logistique ordinale
model_clean <- polr(as.factor(satisfaction) ~ distance_abs_2, data = my_data_frame_clean, method = "logistic")

# Résumé du modèle
summary(model_clean)
Call:
polr(formula = as.factor(satisfaction) ~ distance_abs_2, data = my_data_frame_clean, 
    method = "logistic")

Coefficients:
                  Value Std. Error t value
distance_abs_2 -0.09192    0.03231  -2.845

Intercepts:
            Value    Std. Error t value 
Low|Medium   -0.5649   0.0338   -16.7203
Medium|High   0.5994   0.0339    17.7059

Residual Deviance: 20127.93 
AIC: 20133.93 
code R
# 6. Calculer les probabilités pour chaque niveau de satisfaction
pred_prob_clean <- predict(model_clean, type = "probs")

# Ajouter les probabilités au dataframe nettoyé
my_data_frame_clean$prob_high <- pred_prob_clean[, "High"]
my_data_frame_clean$prob_medium <- pred_prob_clean[, "Medium"]
my_data_frame_clean$prob_low <- pred_prob_clean[, "Low"]

# 7. Visualiser les résultats
ggplot(my_data_frame_clean, aes(x = distance_abs_2)) +
  geom_line(aes(y = prob_high, color = "High"), size = 1) +
  geom_line(aes(y = prob_medium, color = "Medium"), size = 1) +
  geom_line(aes(y = prob_low, color = "Low"), size = 1) +
  labs(title = "Probabilité d'être satisfait selon l'écart aux normes",
       x = "Distance aux Normes (Z-score)",
       y = "Probabilité") +
  scale_color_manual(values = c("High" = "blue", "Medium" = "orange", "Low" = "red")) +
  theme_minimal()

code R
# 1. Charger les librairies nécessaires
library(MASS)
library(ggplot2)

# 2. Vérifier les valeurs manquantes dans le dataframe
#summary(my_data_frame)

# 3. Supprimer les lignes contenant des valeurs manquantes dans les colonnes pertinentes
my_data_frame_clean <- my_data_frame[complete.cases(my_data_frame[c("satisfaction", "Sex", "distance_abs_2")]), ]

# Vérification de la taille du dataframe nettoyé
nrow(my_data_frame_clean)  # Assure-toi que la taille est correcte
[1] 9225
code R
# 4. Convertir les variables catégoriques en facteurs avec les niveaux appropriés
my_data_frame_clean$satisfaction <- factor(my_data_frame_clean$satisfaction, levels = c("Low", "Medium", "High"))
my_data_frame_clean$Sex <- factor(my_data_frame_clean$Sex, levels = c("Men", "Women"))

# 5. Ajuster le modèle de régression logistique ordinale avec interaction entre Sex et distance_abs_2
model_clean_sex <- polr(as.factor(satisfaction) ~ distance_abs_2 * Sex, data = my_data_frame_clean, method = "logistic")

# Résumé du modèle
summary(model_clean_sex)
Call:
polr(formula = as.factor(satisfaction) ~ distance_abs_2 * Sex, 
    data = my_data_frame_clean, method = "logistic")

Coefficients:
                           Value Std. Error t value
distance_abs_2          -0.11679    0.04617 -2.5298
SexWomen                -0.18175    0.06454 -2.8160
distance_abs_2:SexWomen  0.05407    0.06464  0.8365

Intercepts:
            Value    Std. Error t value 
Low|Medium   -0.6607   0.0474   -13.9416
Medium|High   0.5049   0.0472    10.7040

Residual Deviance: 20114.33 
AIC: 20124.33 
code R
# 6. Calculer les probabilités pour chaque niveau de satisfaction par sexe
pred_prob_clean_sex <- predict(model_clean_sex, type = "probs")

# Ajouter les probabilités au dataframe nettoyé
my_data_frame_clean$prob_high <- pred_prob_clean_sex[, "High"]
my_data_frame_clean$prob_medium <- pred_prob_clean_sex[, "Medium"]
my_data_frame_clean$prob_low <- pred_prob_clean_sex[, "Low"]

# 7. Visualiser les résultats en fonction de l'écart aux normes et du sexe (avec facet)
ggplot(my_data_frame_clean, aes(x = distance_abs_2)) +
  geom_line(aes(y = prob_high, color = "High"), size = 1) +
  geom_line(aes(y = prob_medium, color = "Medium"), size = 1) +
  geom_line(aes(y = prob_low, color = "Low"), size = 1) +
  labs(title = "Probabilité d'être satisfait selon l'écart aux normes et par sexe",
       x = "Distance aux Normes (Z-score)",
       y = "Probabilité") +
  scale_color_manual(values = c("High" = "blue", "Medium" = "orange", "Low" = "red")) +
  theme_minimal() +
  facet_wrap(~ Sex) +  # Facette par Sexe
  theme(legend.title = element_blank())  # Enlever le titre de la légende

code R
# 1. Charger les librairies nécessaires
library(MASS)
library(ggplot2)
library(plotly)

# 2. Vérifier les valeurs manquantes dans le dataframe
#summary(my_data_frame)

# 3. Supprimer les lignes contenant des valeurs manquantes dans les colonnes pertinentes
my_data_frame_clean <- my_data_frame[complete.cases(my_data_frame[c("satisfaction", "Sex", "distance_abs_2")]), ]

# Vérification de la taille du dataframe nettoyé
nrow(my_data_frame_clean)  # Assure-toi que la taille est correcte
[1] 9225
code R
# 4. Convertir les variables catégoriques en facteurs avec les niveaux appropriés
my_data_frame_clean$satisfaction <- factor(my_data_frame_clean$satisfaction, levels = c("Low", "Medium", "High"))
my_data_frame_clean$Sex <- factor(my_data_frame_clean$Sex, levels = c("Men", "Women"))

# 5. Ajuster le modèle de régression logistique ordinale avec interaction entre Sex et distance_abs_2
model_clean_sex <- polr(as.factor(satisfaction) ~ distance_abs_2 * Sex, data = my_data_frame_clean, method = "logistic")

# Résumé du modèle
summary(model_clean_sex)
Call:
polr(formula = as.factor(satisfaction) ~ distance_abs_2 * Sex, 
    data = my_data_frame_clean, method = "logistic")

Coefficients:
                           Value Std. Error t value
distance_abs_2          -0.11679    0.04617 -2.5298
SexWomen                -0.18175    0.06454 -2.8160
distance_abs_2:SexWomen  0.05407    0.06464  0.8365

Intercepts:
            Value    Std. Error t value 
Low|Medium   -0.6607   0.0474   -13.9416
Medium|High   0.5049   0.0472    10.7040

Residual Deviance: 20114.33 
AIC: 20124.33 
code R
# 6. Calculer les probabilités pour chaque niveau de satisfaction par sexe
pred_prob_clean_sex <- predict(model_clean_sex, type = "probs")

# Ajouter les probabilités au dataframe nettoyé
my_data_frame_clean$prob_high <- pred_prob_clean_sex[, "High"]
my_data_frame_clean$prob_medium <- pred_prob_clean_sex[, "Medium"]
my_data_frame_clean$prob_low <- pred_prob_clean_sex[, "Low"]

# 7. Visualisation avec ggplot et ajout de commentaires pour les courbes
p <- ggplot(my_data_frame_clean, aes(x = distance_abs_2)) +
  geom_line(aes(y = prob_high, color = "High"), size = 1) +
  geom_line(aes(y = prob_medium, color = "Medium"), size = 1) +
  geom_line(aes(y = prob_low, color = "Low"), size = 1) +
  labs(title = "Probabilité d'être satisfait selon l'écart aux normes et par sexe",
       x = "Distance aux Normes (Z-score)",
       y = "Probabilité") +
  scale_color_manual(values = c("High" = "blue", "Medium" = "orange", "Low" = "red")) +
  theme_minimal() +
  facet_wrap(~ Sex) +  # Facette par Sexe
  theme(legend.title = element_blank())  # Enlever le titre de la légende

# 8. Rendre le graphique interactif avec plotly et ajouter des annotations
p_interactive <- ggplotly(p)

# Ajouter des annotations pour chaque courbe
p_interactive <- p_interactive %>%
  layout(
    annotations = list(
      list(
        x = 0.5, y = 0.8, 
        text = "Probabilité High", 
        showarrow = TRUE, 
        arrowhead = 2,
        ax = -50, ay = -50
      ),
      list(
        x = 0.5, y = 0.6, 
        text = "Probabilité Medium", 
        showarrow = TRUE, 
        arrowhead = 2,
        ax = -50, ay = -50
      ),
      list(
        x = 0.5, y = 0.4, 
        text = "Probabilité Low", 
        showarrow = TRUE, 
        arrowhead = 2,
        ax = -50, ay = -50
      )
    )
  )

# Afficher le graphique interactif
p_interactive
code R
# 1. Charger les librairies nécessaires
library(MASS)
library(ggplot2)
library(plotly)

# 2. Convertir la colonne de réponses A15 en une variable Santé catégorique (3 niveaux)
# Supposons que ta variable santé est dans la colonne 'A15' de 'my_data_frame'

# Remplacer les valeurs de A15 par les niveaux correspondants pour Santé
my_data_frame$Santé <- factor(my_data_frame$A15,
                              levels = c(1, 2, 3, 4, 5),
                              labels = c("Very Good", "Good", "Fair", "Poor", "Very Poor"))

# Convertir en 3 niveaux: "Good", "Fair", "Poor"
my_data_frame$Santé <- factor(my_data_frame$Santé,
                              levels = c("Very Good", "Good", "Fair", "Poor", "Very Poor"),
                              labels = c("Good", "Good", "Fair", "Poor", "Poor"))

# Exclure les réponses "Ne sait pas" (6) et "Refus" (7)
my_data_frame_clean <- my_data_frame[!(my_data_frame$A15 %in% c(6, 7)), ]

# 3. Vérifier que la variable Santé est correctement créée
table(my_data_frame_clean$Santé)

Good Fair Poor 
6410 1998  783 
code R
# 4. Convertir 'Sex' en facteur
my_data_frame_clean$Sex <- factor(my_data_frame_clean$Sex, levels = c("Men", "Women"))

# 5. Ajuster le modèle de régression logistique ordinale avec interaction entre Sex et distance_abs_2
model_clean_sex_health <- polr(Santé ~ distance_abs_2 * Sex, data = my_data_frame_clean, method = "logistic")

# Résumé du modèle
summary(model_clean_sex_health)
Call:
polr(formula = Santé ~ distance_abs_2 * Sex, data = my_data_frame_clean, 
    method = "logistic")

Coefficients:
                           Value Std. Error t value
distance_abs_2          -0.19064    0.05772 -3.3026
SexWomen                 0.02624    0.07578  0.3463
distance_abs_2:SexWomen  0.21223    0.07738  2.7426

Intercepts:
          Value   Std. Error t value
Good|Fair  0.7948  0.0554    14.3528
Fair|Poor  2.3365  0.0628    37.2049

Residual Deviance: 14545.71 
AIC: 14555.71 
code R
# 6. Calculer les probabilités pour chaque niveau de Santé par sexe
pred_prob_clean_sex_health <- predict(model_clean_sex_health, type = "probs")

# Ajouter les probabilités au dataframe nettoyé
my_data_frame_clean$prob_poor <- pred_prob_clean_sex_health[, "Poor"]
my_data_frame_clean$prob_fair <- pred_prob_clean_sex_health[, "Fair"]
my_data_frame_clean$prob_good <- pred_prob_clean_sex_health[, "Good"]

# 7. Visualisation avec ggplot et ajout de commentaires pour les courbes
p_health <- ggplot(my_data_frame_clean, aes(x = distance_abs_2)) +
  geom_line(aes(y = prob_poor, color = "Poor"), size = 1) +
  geom_line(aes(y = prob_fair, color = "Fair"), size = 1) +
  geom_line(aes(y = prob_good, color = "Good"), size = 1) +
  labs(title = "Probabilité d'être en bonne santé selon l'écart aux normes et par sexe",
       x = "Distance aux Normes (Z-score)",
       y = "Probabilité") +
  scale_color_manual(values = c("Poor" = "red", 
                                "Fair" = "orange", 
                                "Good" = "green")) +
  theme_minimal() +
  facet_wrap(~ Sex) +  # Facette par Sexe
  theme(legend.title = element_blank())  # Enlever le titre de la légende

# 8. Rendre le graphique interactif avec plotly et ajouter des annotations
p_interactive_health <- ggplotly(p_health)

# Ajouter des annotations pour chaque courbe
p_interactive_health <- p_interactive_health %>%
  layout(
    annotations = list(
      list(
        x = 0.5, y = 0.8, 
        text = "Probabilité Poor", 
        showarrow = TRUE, 
        arrowhead = 2,
        ax = -50, ay = -50
      ),
      list(
        x = 0.5, y = 0.6, 
        text = "Probabilité Fair", 
        showarrow = TRUE, 
        arrowhead = 2,
        ax = -50, ay = -50
      ),
      list(
        x = 0.5, y = 0.4, 
        text = "Probabilité Good", 
        showarrow = TRUE, 
        arrowhead = 2,
        ax = -50, ay = -50
      )
    )
  )

# Afficher le graphique interactif
p_interactive_health
code R
# Créer la variable revenu en fonction de CRITREVENU
my_data_frame$revenu <- NA  # Initialisation de la variable

# Assignation des tranches de revenu
my_data_frame$revenu[my_data_frame$CRITREVENU %in% c(1, 2)] <- "Moins de 1000 euros"
my_data_frame$revenu[my_data_frame$CRITREVENU %in% c(3, 4)] <- "1000 à 1499 euros"
my_data_frame$revenu[my_data_frame$CRITREVENU %in% c(5, 6)] <- "1500 à 2499 euros"
my_data_frame$revenu[my_data_frame$CRITREVENU %in% c(7, 8, 9, 10)] <- "2500 euros et plus"

# Convertir la variable revenu en facteur
my_data_frame$revenu <- factor(my_data_frame$revenu, 
                               levels = c("Moins de 1000 euros", "1000 à 1499 euros", "1500 à 2499 euros", "2500 euros et plus"))

# Vérifier les niveaux de la nouvelle variable revenu
summary(my_data_frame$revenu)
Moins de 1000 euros   1000 à 1499 euros   1500 à 2499 euros  2500 euros et plus 
                766                1280                2241                3775 
               NA's 
               1172 
code R
# Charger les librairies nécessaires
library(MASS)
library(ggplot2)
library(plotly)

# 2. Vérifier les valeurs manquantes dans le dataframe
#summary(my_data_frame)

# 3. Supprimer les lignes contenant des valeurs manquantes dans les colonnes pertinentes
my_data_frame_clean <- my_data_frame[complete.cases(my_data_frame[c("satisfaction", "Sex", "distance_abs_2", "revenu")]), ]

# Vérification de la taille du dataframe nettoyé
nrow(my_data_frame_clean)  # Assure-toi que la taille est correcte
[1] 8057
code R
# 4. Convertir les variables catégoriques en facteurs avec les niveaux appropriés
my_data_frame_clean$satisfaction <- factor(my_data_frame_clean$satisfaction, levels = c("Low", "Medium", "High"))
my_data_frame_clean$Sex <- factor(my_data_frame_clean$Sex, levels = c("Men", "Women"))


my_data_frame_clean$revenu <- factor(my_data_frame_clean$revenu, 
                                     levels = c("Moins de 1000 euros", "1000 à 1499 euros", 
                                                "1500 à 2499 euros", "2500 euros et plus"),
                                     labels = c("<1000", "[1000-1499]", "[1500-2499]", "≥2500"))

# 5. Ajuster le modèle de régression logistique ordinale avec interaction entre Sexe, revenu et distance_abs_2
model_clean_sex_revenu <- polr(as.factor(satisfaction) ~ distance_abs_2 * Sex * revenu, 
                               data = my_data_frame_clean, method = "logistic")

# Résumé du modèle
summary(model_clean_sex_revenu)
Call:
polr(formula = as.factor(satisfaction) ~ distance_abs_2 * Sex * 
    revenu, data = my_data_frame_clean, method = "logistic")

Coefficients:
                                             Value Std. Error t value
distance_abs_2                            -0.15950     0.1647 -0.9683
SexWomen                                  -0.32601     0.2353 -1.3856
revenu[1000-1499]                         -0.06179     0.2271 -0.2721
revenu[1500-2499]                         -0.42024     0.1997 -2.1042
revenu≥2500                               -1.08504     0.1889 -5.7434
distance_abs_2:SexWomen                    0.22862     0.2295  0.9961
distance_abs_2:revenu[1000-1499]          -0.02571     0.2198 -0.1170
distance_abs_2:revenu[1500-2499]          -0.11661     0.1912 -0.6099
distance_abs_2:revenu≥2500                 0.14320     0.1780  0.8046
SexWomen:revenu[1000-1499]                 0.04886     0.2999  0.1629
SexWomen:revenu[1500-2499]                -0.06160     0.2701 -0.2280
SexWomen:revenu≥2500                       0.07601     0.2558  0.2971
distance_abs_2:SexWomen:revenu[1000-1499] -0.09139     0.2959 -0.3089
distance_abs_2:SexWomen:revenu[1500-2499]  0.03572     0.2659  0.1343
distance_abs_2:SexWomen:revenu≥2500       -0.32040     0.2501 -1.2813

Intercepts:
            Value   Std. Error t value
Low|Medium  -1.3544  0.1766    -7.6705
Medium|High -0.1353  0.1759    -0.7691

Residual Deviance: 17140.82 
AIC: 17174.82 
code R
# 6. Calculer les probabilités pour chaque niveau de satisfaction par sexe et revenu
pred_prob_clean_sex_revenu <- predict(model_clean_sex_revenu, type = "probs")

# Ajouter les probabilités au dataframe nettoyé
my_data_frame_clean$prob_high <- pred_prob_clean_sex_revenu[, "High"]
my_data_frame_clean$prob_medium <- pred_prob_clean_sex_revenu[, "Medium"]
my_data_frame_clean$prob_low <- pred_prob_clean_sex_revenu[, "Low"]

p <- ggplot(my_data_frame_clean, aes(x = distance_abs_2)) +
  geom_line(aes(y = prob_high, color = "High"), size = 1) +
  geom_line(aes(y = prob_medium, color = "Medium"), size = 1) +
  geom_line(aes(y = prob_low, color = "Low"), size = 1) +
  labs(title = "Probabilité d'être satisfait selon l'écart aux normes, par Sexe et Revenu",
       x = "Distance aux Normes (Z-score)",
       y = "Probabilité") +
  scale_color_manual(values = c("High" = "blue", "Medium" = "orange", "Low" = "red")) +
  theme_minimal() +
  facet_grid(revenu ~ Sex) +  # Facet verticalement par revenu et horizontalement par sexe
  theme(
    legend.title = element_blank(),
    strip.text.y = element_text(size = 10),  # Réduire la taille des labels des facettes
    panel.spacing = unit(2, "lines")  # Augmenter l'espacement vertical entre facettes
  )

# Rendre interactif
p_interactive <- ggplotly(p)

# Afficher le graphique
p_interactive
code R
# Charger les librairies nécessaires
library(MASS)
library(ggplot2)
library(plotly)



# Supprimer les lignes contenant des valeurs manquantes dans les colonnes pertinentes
my_data_frame_clean <- my_data_frame[complete.cases(my_data_frame[c("Santé", "Sex", "distance_abs_2", "revenu")]), ]

# Convertir les variables catégoriques en facteurs avec les niveaux appropriés
my_data_frame_clean$Sex <- factor(my_data_frame_clean$Sex, levels = c("Men", "Women"))

my_data_frame_clean$revenu <- factor(my_data_frame_clean$revenu, 
                                     levels = c("Moins de 1000 euros", "1000 à 1499 euros", 
                                                "1500 à 2499 euros", "2500 euros et plus"),
                                     labels = c("<1000", "[1000-1499]", "[1500-2499]", "≥2500"))

# Ajuster le modèle de régression logistique ordinale avec interaction entre Sexe, revenu et distance_abs_2
model_clean_sex_revenu <- polr(Santé ~ distance_abs_2 * Sex * revenu, 
                               data = my_data_frame_clean, method = "logistic")

# Résumé du modèle
summary(model_clean_sex_revenu)
Call:
polr(formula = Santé ~ distance_abs_2 * Sex * revenu, data = my_data_frame_clean, 
    method = "logistic")

Coefficients:
                                             Value Std. Error t value
distance_abs_2                            -0.07229     0.1772 -0.4079
SexWomen                                   0.02646     0.2401  0.1102
revenu[1000-1499]                         -0.17956     0.2341 -0.7672
revenu[1500-2499]                         -0.61613     0.2096 -2.9391
revenu≥2500                               -1.35968     0.2034 -6.6848
distance_abs_2:SexWomen                    0.09937     0.2390  0.4157
distance_abs_2:revenu[1000-1499]          -0.08366     0.2352 -0.3558
distance_abs_2:revenu[1500-2499]          -0.20282     0.2118 -0.9578
distance_abs_2:revenu≥2500                -0.11564     0.2024 -0.5713
SexWomen:revenu[1000-1499]                 0.09874     0.3086  0.3199
SexWomen:revenu[1500-2499]                -0.29406     0.2841 -1.0349
SexWomen:revenu≥2500                      -0.07925     0.2746 -0.2886
distance_abs_2:SexWomen:revenu[1000-1499] -0.11806     0.3100 -0.3809
distance_abs_2:SexWomen:revenu[1500-2499]  0.26781     0.2862  0.9357
distance_abs_2:SexWomen:revenu≥2500        0.15741     0.2737  0.5750

Intercepts:
          Value   Std. Error t value
Good|Fair -0.0474  0.1802    -0.2629
Fair|Poor  1.5620  0.1821     8.5795

Residual Deviance: 12216.75 
AIC: 12250.75 
code R
# Calculer les probabilités pour chaque niveau de santé par sexe et revenu
pred_prob_clean_sex_revenu <- predict(model_clean_sex_revenu, type = "probs")

# Ajouter les probabilités au dataframe nettoyé
my_data_frame_clean$prob_high <- pred_prob_clean_sex_revenu[, "Good"]
my_data_frame_clean$prob_medium <- pred_prob_clean_sex_revenu[, "Fair"]
my_data_frame_clean$prob_low <- pred_prob_clean_sex_revenu[, "Poor"]

# Visualisation avec ggplot
p <- ggplot(my_data_frame_clean, aes(x = distance_abs_2)) +
  geom_line(aes(y = prob_high, color = "High"), size = 1) +
  geom_line(aes(y = prob_medium, color = "Medium"), size = 1) +
  geom_line(aes(y = prob_low, color = "Low"), size = 1) +
  labs(title = "Probabilité d'avoir un bon état de santé selon l'écart aux normes, par Sexe et Revenu",
       x = "Distance aux Normes (Z-score)",
       y = "Probabilité") +
  scale_color_manual(values = c("High" = "blue", "Medium" = "orange", "Low" = "red")) +
  theme_minimal() +
  facet_grid(revenu ~ Sex) +  # Facet verticalement par revenu et horizontalement par sexe
  theme(
    legend.title = element_blank(),
    strip.text.y = element_text(size = 10),  # Réduire la taille des labels des facettes
    panel.spacing = unit(2, "lines")  # Augmenter l'espacement vertical entre facettes
  )

# Rendre interactif
p_interactive <- ggplotly(p)

# Afficher le graphique
p_interactive
code R
my_data_frame$distance_abs_2_bins<-as.factor(my_data_frame$distance_abs_2_bins)
my_data_frame$satisfaction<-as.factor(my_data_frame$satisfaction)
rego <- MASS::polr(
  satisfaction ~ Sex + revenu + age_group + distance_abs_2_bins,
  data = my_data_frame
)

theme_gtsummary_language("en", decimal.mark = ",")
rego |> 
  tbl_regression(
    exponentiate = TRUE,
    tidy_fun = broom.helpers::tidy_parameters
  ) |> 
  bold_labels() |> 
  add_global_p(keep = TRUE) |> 
  as_kable_extra(format = "latex", booktabs = TRUE) |> 
  kable_styling(latex_options = c("hold_position", "scale_down"), full_width = FALSE)
code R
rego |> 
  ggstats::ggcoef_table(
    exponentiate = TRUE,
    tidy_fun = broom.helpers::tidy_parameters
  )

code R
rego |> 
  broom.helpers::plot_marginal_predictions() |> 
  patchwork::wrap_plots(ncol = 1) &
  scale_y_continuous(labels = scales::percent) &
  coord_flip()

code R
library(dplyr)

# Vérifier et nettoyer la colonne Sex
my_data_frame <- my_data_frame %>%
  mutate(
    Sex = trimws(as.character(Sex)),  # Supprime les espaces autour des valeurs
    Sex = case_when(
      tolower(Sex) == "men" ~ "Men",   # Uniformise "men" en "Men"
      tolower(Sex) == "women" ~ "Women",  # Uniformise "women" en "Women"
      TRUE ~ NA_character_  # Remplace toutes les autres valeurs par NA
    )
  ) %>%
  drop_na(Sex)  # Supprime les lignes où Sex est NA

# Séparer les données
df_men <- my_data_frame %>% filter(Sex == "Men")
df_women <- my_data_frame %>% filter(Sex == "Women")

# Vérifier la séparation
print(dim(df_men))   # Nombre de lignes et colonnes pour les hommes
[1] 4162 1593
code R
print(dim(df_women)) # Nombre de lignes et colonnes pour les femmes
[1] 5072 1593
code R
# Modèles de régression
rego_men <- polr(
  satisfaction ~ revenu + age_group + score_scale,
  data = df_men
)



# Fonction pour générer le tableau
generate_table <- function(rego_model) {
  rego_model |> 
    tbl_regression(
      exponentiate = TRUE,
      tidy_fun = broom.helpers::tidy_parameters
    ) |> 
    bold_labels() |> 
    add_global_p(keep = TRUE) |> 
    as_kable_extra(format = "latex", booktabs = TRUE) |> 
    kable_styling(latex_options = c("hold_position", "scale_down"), full_width = FALSE)
}

p1<- rego_men |> 
  ggstats::ggcoef_table(
    exponentiate = TRUE,
    tidy_fun = broom.helpers::tidy_parameters
  )

p1

code R
p2<-rego_men |> 
  broom.helpers::plot_marginal_predictions() |> 
  patchwork::wrap_plots(ncol = 1) &
  scale_y_continuous(labels = scales::percent) &
  coord_flip()
p2

code R
df_women <- df_women %>%
  mutate(across(where(is.factor), droplevels))
rego_women <- polr(
  satisfaction ~ revenu + age_group + score_scale,
  data = df_women
)

summary(rego_women)
Call:
polr(formula = satisfaction ~ revenu + age_group + score_scale, 
    data = df_women)

Coefficients:
                           Value Std. Error  t value
revenu1000 à 1499 euros   0.2604    0.11542   2.2566
revenu1500 à 2499 euros   0.4406    0.10790   4.0832
revenu2500 euros et plus  0.5379    0.10226   5.2603
age_group[38-54[         -0.2421    0.07518  -3.2199
age_group[54-67[         -0.8566    0.07840 -10.9250
age_group[67-97[         -1.8190    0.11709 -15.5350
score_scale2             -0.2358    0.47205  -0.4995
score_scale3             -0.2969    0.46379  -0.6402
score_scale4             -0.1785    0.46406  -0.3847
score_scale5             -0.2060    0.46711  -0.4410
score_scaleVery Feminine -0.1623    0.49138  -0.3302

Intercepts:
           Value    Std. Error t value 
High|Low    -1.1607   0.4724    -2.4567
Low|Medium   0.6977   0.4721     1.4778

Residual Deviance: 9113.95 
AIC: 9139.95 
(682 observations effacées parce que manquantes)
code R
p3<-rego_women |> 
  ggstats::ggcoef_table(
    exponentiate = TRUE,
    tidy_fun = broom.helpers::tidy_parameters
  )
p3

code R
p4<-rego_women |> 
  broom.helpers::plot_marginal_predictions() |> 
  patchwork::wrap_plots(ncol = 1) &
  scale_y_continuous(labels = scales::percent) &
  coord_flip()

p4

code R
library(dplyr)
library(ggplot2)
library(plotly)

# Créer un tableau de proportions
df_proportions <- my_data_frame %>%
  group_by(revenu, score_scale) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(revenu) %>%
  mutate(proportion = count / sum(count))  # Calcul de la proportion

# Graphique interactif avec Plotly (barres empilées)
fig <- plot_ly(df_proportions, 
               x = ~revenu, 
               y = ~proportion, 
               color = ~score_scale, 
               type = "bar",
               text = ~paste0(round(proportion*100, 1), "%"), 
               textposition = "inside") %>%
  layout(title = "Proportion de score_scale par tranche de revenu",
         xaxis = list(title = "Tranche de revenu"),
         yaxis = list(title = "Proportion", tickformat = "%"),
         barmode = "stack")  # Empilement des barres

fig  # Affichage du graphique interactif
code R
library(dplyr)
library(ggplot2)
library(plotly)

# Créer un tableau de proportions
df_proportions <- my_data_frame %>%
  group_by(CLASSIF, score_scale) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(CLASSIF) %>%
  mutate(proportion = count / sum(count))  # Calcul de la proportion

# Graphique interactif avec Plotly (barres empilées)
fig <- plot_ly(df_proportions, 
               x = ~CLASSIF, 
               y = ~proportion, 
               color = ~score_scale, 
               type = "bar",
               text = ~paste0(round(proportion*100, 1), "%"), 
               textposition = "inside") %>%
  layout(title = "Proportion de score_scale par profession",
         xaxis = list(title = "Profession"),
         yaxis = list(title = "Proportion", tickformat = "%"),
         barmode = "stack")  # Empilement des barres

fig  # Affichage du graphique interactif
1 Manoeuvre ou ouvrier spécialisé
2 Ouvrier qualifié ou hautement qualifié/ technicien(ne) d’atelier
3 Technicien(ne)
4 Agent de maîtrise, maîtrise administrative ou commerciale, VRP (non cadre)
5 Ingénieur, Cadre
6 Employé(e) de bureau, Employé(e) de commerce, Personnel de services
7 Directeur général, Adjoint direct
8 NSP
9 REF
code R
library(dplyr)
library(ggplot2)
library(plotly)

# Créer un tableau de proportions
df_proportions <- my_data_frame %>%
  group_by(satisfaction, score_scale) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(satisfaction) %>%
  mutate(proportion = count / sum(count))  # Calcul de la proportion

# Graphique interactif avec Plotly (barres empilées)
fig <- plot_ly(df_proportions, 
               x = ~satisfaction, 
               y = ~proportion, 
               color = ~score_scale, 
               type = "bar",
               text = ~paste0(round(proportion*100, 1), "%"), 
               textposition = "inside") %>%
  layout(title = "Proportion de score_scale par degré de satisfaction",
         xaxis = list(title = "Satisfaction"),
         yaxis = list(title = "Proportion", tickformat = "%"),
         barmode = "stack")  # Empilement des barres

fig  # Affichage du graphique interactif
code R
library(dplyr)
library(ggplot2)
library(plotly)

# Convertir DIPLOM en facteur
my_data_frame$DIPLOM <- as.factor(my_data_frame$DIPLOM)

# Créer un tableau de proportions
df_proportions <- my_data_frame %>%
  group_by(DIPLOM, score_scale) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(DIPLOM) %>%
  mutate(proportion = count / sum(count))  # Calcul de la proportion

# Ajouter une colonne avec les descriptions des diplômes
diplome_labels <- c(
  "Vous n'avez jamais été à l'école ou vous l'avez quittée avant la fin du primaire",
  "Aucun diplôme et scolarité interrompue à la fin du primaire ou avant la fin du collège",
  "Aucun diplôme et scolarité jusqu'à la fin du collège et au-delà",
  "CEP",
  "BEPC, brevet élémentaire, brevet des collèges, DNB",
  "CAP, BEP ou diplôme équivalent",
  "Baccalauréat général ou technologique, brevet supérieur",
  "Capacité en droit, DAEU, ESEU",
  "Baccalauréat professionnel, brevet professionnel, de technicien ou d'enseignement, diplôme équivalent",
  "BTS, DUT, DEUST, diplôme de la santé ou social de niveau Bac+2 ou diplôme équivalent",
  "Licence, licence pro, maîtrise ou autre diplôme de niveau Bac+3 ou 4 ou diplôme équivalent",
  "Master, DEA, DESS, diplôme grande école de niveau Bac+5, doctorat de santé",
  "Doctorat de recherche (hors santé)",
  "NSP",
  "REF"
)

# Ajouter les libellés des diplômes à df_proportions
df_proportions <- df_proportions %>%
  mutate(DIPLOM_label = diplome_labels[as.numeric(DIPLOM)])

# Graphique interactif avec Plotly (barres empilées)
fig <- plot_ly(df_proportions, 
               x = ~DIPLOM, 
               y = ~proportion, 
               color = ~score_scale, 
               type = "bar",
               text = ~paste(DIPLOM_label, "<br>", round(proportion*100, 1), "%"),  # Affichage du libellé et proportion
               textposition = "inside") %>%
  layout(title = "Proportion de score_scale par niveau de diplôme",
         xaxis = list(title = "Niveau de diplôme"),
         yaxis = list(title = "Proportion", tickformat = "%"),
         barmode = "stack",  # Empilement des barres
         hovermode = "closest")  # Afficher les informations les plus proches du survol

fig
code R
library(dplyr)
library(ggplot2)
library(plotly)

# Convertir DIPLOM en facteur
my_data_frame$DIPLOM <- as.factor(my_data_frame$DIPLOM)
# Création des quartiles de distance_abs_2
my_data_frame <- my_data_frame %>%
  mutate(quartile_distance = cut(distance_abs_2, 
                                 breaks = quantile(distance_abs_2, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE), 
                                 include.lowest = TRUE, 
                                 labels = c("Q1", "Q2", "Q3", "Q4")))  # Étiquettes des quartiles

# Créer un tableau de proportions
df_proportions <- my_data_frame %>%
  group_by(DIPLOM, quartile_distance) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(DIPLOM) %>%
  mutate(proportion = count / sum(count))  # Calcul de la proportion

# Ajouter une colonne avec les descriptions des diplômes
diplome_labels <- c(
  "Vous n'avez jamais été à l'école ou vous l'avez quittée avant la fin du primaire",
  "Aucun diplôme et scolarité interrompue à la fin du primaire ou avant la fin du collège",
  "Aucun diplôme et scolarité jusqu'à la fin du collège et au-delà",
  "CEP",
  "BEPC, brevet élémentaire, brevet des collèges, DNB",
  "CAP, BEP ou diplôme équivalent",
  "Baccalauréat général ou technologique, brevet supérieur",
  "Capacité en droit, DAEU, ESEU",
  "Baccalauréat professionnel, brevet professionnel, de technicien ou d'enseignement, diplôme équivalent",
  "BTS, DUT, DEUST, diplôme de la santé ou social de niveau Bac+2 ou diplôme équivalent",
  "Licence, licence pro, maîtrise ou autre diplôme de niveau Bac+3 ou 4 ou diplôme équivalent",
  "Master, DEA, DESS, diplôme grande école de niveau Bac+5, doctorat de santé",
  "Doctorat de recherche (hors santé)",
  "NSP",
  "REF"
)

# Ajouter les libellés des diplômes à df_proportions
df_proportions <- df_proportions %>%
  mutate(DIPLOM_label = diplome_labels[as.numeric(DIPLOM)])

# Graphique interactif avec Plotly (barres empilées)
fig <- plot_ly(df_proportions, 
               x = ~DIPLOM, 
               y = ~proportion, 
               color = ~quartile_distance, 
               type = "bar",
               text = ~paste(DIPLOM_label, "<br>", round(proportion*100, 1), "%"),  # Affichage du libellé et proportion
               textposition = "inside") %>%
  layout(title = "Quartile Distance par niveau de diplôme",
         xaxis = list(title = "Niveau de diplôme"),
         yaxis = list(title = "Proportion", tickformat = "%"),
         barmode = "stack",  # Empilement des barres
         hovermode = "closest")  # Afficher les informations les plus proches du survol

fig
code R
library(dplyr)
library(ggplot2)
library(plotly)

# Création des quartiles de distance_abs_2
my_data_frame <- my_data_frame %>%
  mutate(quartile_distance = cut(distance_abs_2, 
                                 breaks = quantile(distance_abs_2, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE), 
                                 include.lowest = TRUE, 
                                 labels = c("Q1", "Q2", "Q3", "Q4")))  # Étiquettes des quartiles

table(my_data_frame$quartile_distance)  # Vérification

  Q1   Q2   Q3   Q4 
2309 2308 2308 2309 
code R
# Calcul des proportions correctes par satisfaction (A2)
df_proportions <- my_data_frame %>%
  group_by(satisfaction, quartile_distance) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(satisfaction) %>%  # Normalisation par A2 et non par quartile_distance
  mutate(proportion = count / sum(count))  

# Graphique interactif avec Plotly (barres empilées)
fig <- plot_ly(df_proportions, 
               x = ~satisfaction, 
               y = ~proportion, 
               color = ~quartile_distance, 
               type = "bar",
               text = ~paste0(round(proportion * 100, 1), "%"), 
               textposition = "inside") %>%
  layout(title = "Proportion de distance aux normes par degré de satisfaction",
         xaxis = list(title = "Satisfaction"),
         yaxis = list(title = "Proportion", tickformat = "%"),
         barmode = "stack")  # Empilement des barres pour que chaque A2 fasse 100%

fig
code R
# Créer un tableau de contingence
contingency_table <- table(my_data_frame$satisfaction, my_data_frame$quartile_distance)

# Test du chi-deux
chi2_test <- chisq.test(contingency_table)

# Résultats
chi2_test

    Pearson's Chi-squared test

data:  contingency_table
X-squared = 11.811, df = 6, p-value = 0.06633
code R
contingency_table <- table(my_data_frame$satisfaction, my_data_frame$quartile_distance)
print(contingency_table)
        
          Q1  Q2  Q3  Q4
  High   785 779 823 730
  Low    848 869 861 925
  Medium 674 657 621 653
code R
library(nnet)

my_data_frame$satisfaction <- as.factor(my_data_frame$satisfaction)  # S'assurer que A2 est factorielle
my_data_frame$quartile_distance <- as.factor(my_data_frame$quartile_distance)  # Idem

model <- multinom(satisfaction ~ quartile_distance, data = my_data_frame)
# weights:  15 (8 variable)
initial  value 10134.698363 
iter  10 value 10062.232209
final  value 10062.112229 
converged
code R
summary(model)
Call:
multinom(formula = satisfaction ~ quartile_distance, data = my_data_frame)

Coefficients:
       (Intercept) quartile_distanceQ2 quartile_distanceQ3 quartile_distanceQ4
Low     0.07719688          0.03213488         -0.03205884          0.15955098
Medium -0.15245629         -0.01787090         -0.12916843          0.04098897

Std. Errors:
       (Intercept) quartile_distanceQ2 quartile_distanceQ3 quartile_distanceQ4
Low     0.04952907          0.06991118          0.06949558          0.07002909
Medium  0.05251254          0.07458776          0.07471907          0.07522514

Residual Deviance: 20124.22 
AIC: 20140.22 
code R
library(ggeffects)
pred <- ggeffect(model, terms = "quartile_distance")
plot(pred)  # Visualisation des probabilités

1.4 Diplomes et distance aux normes de genre

code R
library(dplyr)
library(ggplot2)
library(plotly)

# Convertir DIPLOM et Sex en facteurs
my_data_frame$DIPLOM <- as.factor(my_data_frame$DIPLOM)
my_data_frame$Sex <- as.factor(my_data_frame$Sex)

# Création des quartiles de distance_abs_2
my_data_frame <- my_data_frame %>%
  mutate(quartile_distance = cut(distance_abs_2, 
                                 breaks = quantile(distance_abs_2, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE), 
                                 include.lowest = TRUE, 
                                 labels = c("Q1", "Q2", "Q3", "Q4")))  # Étiquettes des quartiles

# Créer un tableau de proportions
df_proportions <- my_data_frame %>%
  group_by(Sex, DIPLOM, quartile_distance) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(Sex, DIPLOM) %>%
  mutate(proportion = count / sum(count))  # Calcul de la proportion

# Ajouter une colonne avec les descriptions des diplômes
diplome_labels <- c(
  "Vous n'avez jamais été à l'école ou vous l'avez quittée avant la fin du primaire",
  "Aucun diplôme et scolarité interrompue à la fin du primaire ou avant la fin du collège",
  "Aucun diplôme et scolarité jusqu'à la fin du collège et au-delà",
  "CEP",
  "BEPC, brevet élémentaire, brevet des collèges, DNB",
  "CAP, BEP ou diplôme équivalent",
  "Baccalauréat général ou technologique, brevet supérieur",
  "Capacité en droit, DAEU, ESEU",
  "Baccalauréat professionnel, brevet professionnel, de technicien ou d'enseignement, diplôme équivalent",
  "BTS, DUT, DEUST, diplôme de la santé ou social de niveau Bac+2 ou diplôme équivalent",
  "Licence, licence pro, maîtrise ou autre diplôme de niveau Bac+3 ou 4 ou diplôme équivalent",
  "Master, DEA, DESS, diplôme grande école de niveau Bac+5, doctorat de santé",
  "Doctorat de recherche (hors santé)",
  "NSP",
  "REF"
)

# Ajouter les libellés des diplômes à df_proportions
df_proportions <- df_proportions %>%
  mutate(DIPLOM_label = diplome_labels[as.numeric(DIPLOM)])

# Séparer les données en deux sous-ensembles (Hommes et Femmes)
df_men <- df_proportions %>% filter(Sex == "Men")
df_women <- df_proportions %>% filter(Sex == "Women")

# Graphique pour les hommes
fig_men <- plot_ly(df_men, 
                   x = ~DIPLOM, 
                   y = ~proportion, 
                   color = ~quartile_distance, 
                   type = "bar",
                   text = ~paste(DIPLOM_label, "<br>", round(proportion*100, 1), "%"),  # Affichage du libellé et proportion
                   textposition = "inside") %>%
  layout(title = "Proportion de score_scale par niveau de diplôme (Hommes)",
         xaxis = list(title = "Niveau de diplôme"),
         yaxis = list(title = "Proportion", tickformat = "%"),
         barmode = "stack",  # Empilement des barres
         hovermode = "closest")  # Afficher les informations les plus proches du survol

# Graphique pour les femmes
fig_women <- plot_ly(df_women, 
                     x = ~DIPLOM, 
                     y = ~proportion, 
                     color = ~quartile_distance, 
                     type = "bar",
                     text = ~paste(DIPLOM_label, "<br>", round(proportion*100, 1), "%"),  # Affichage du libellé et proportion
                     textposition = "inside") %>%
  layout(title = "Proportion de score_scale par niveau de diplôme (Femmes)",
         xaxis = list(title = "Niveau de diplôme"),
         yaxis = list(title = "Proportion", tickformat = "%"),
         barmode = "stack",  # Empilement des barres
         hovermode = "closest")  # Afficher les informations les plus proches du survol

# Afficher les deux graphiques séparés
fig_men
code R
fig_women
code R
library(dplyr)
library(ggplot2)
library(plotly)

# Création des quartiles de distance_abs_2
my_data_frame <- my_data_frame %>%
  mutate(quartile_distance = cut(distance_abs_2, 
                                 breaks = quantile(distance_abs_2, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE), 
                                 include.lowest = TRUE, 
                                 labels = c("Q1", "Q2", "Q3", "Q4")))  # Étiquettes des quartiles

table(my_data_frame$quartile_distance)  # Vérification

  Q1   Q2   Q3   Q4 
2309 2308 2308 2309 
code R
# Calcul des proportions correctes par satisfaction (A2)
df_proportions <- my_data_frame %>%
  group_by(Santé, quartile_distance) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(Santé) %>%  # Normalisation par A2 et non par quartile_distance
  mutate(proportion = count / sum(count))  

# Graphique interactif avec Plotly (barres empilées)
fig <- plot_ly(df_proportions, 
               x = ~Santé, 
               y = ~proportion, 
               color = ~quartile_distance, 
               type = "bar",
               text = ~paste0(round(proportion * 100, 1), "%"), 
               textposition = "inside") %>%
  layout(title = "Proportion de distance aux normes par état de Santé",
         xaxis = list(title = "Santé"),
         yaxis = list(title = "Proportion", tickformat = "%"),
         barmode = "stack")  # Empilement des barres pour que chaque A2 fasse 100%

fig
code R
# Créer un tableau de contingence
contingency_table <- table(my_data_frame$Santé, my_data_frame$quartile_distance)

# Test du chi-deux
chi2_test <- chisq.test(contingency_table)

1.5 Comparaison des modèles (Suite)

code R
my_data_frame$Sex <- as.factor(my_data_frame$Sex)

# Modèle avec 'identity'
model_identity <- glm(Sex ~ identity, data = my_data_frame, family = binomial)

# Modèle avec 'score'
model_score <- glm(Sex ~ score, data = my_data_frame, family = binomial)

# Comparer les modèles via les critères AIC (Akaike Information Criterion)
aic_identity <- AIC(model_identity)
aic_score <- AIC(model_score)

# Comparer les modèles via la deviance
deviance_identity <- deviance(model_identity)
deviance_score <- deviance(model_score)

# Résumé des modèles
summary(model_identity)

Call:
glm(formula = Sex ~ identity, family = binomial, data = my_data_frame)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.00930    0.08534   35.26   <2e-16 ***
identity    -7.59704    0.22373  -33.96   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12711  on 9233  degrees of freedom
Residual deviance: 11162  on 9232  degrees of freedom
AIC: 11166

Number of Fisher Scoring iterations: 4
code R
summary(model_score)

Call:
glm(formula = Sex ~ score, family = binomial, data = my_data_frame)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.58503    0.03010  -19.43   <2e-16 ***
score        1.04195    0.02254   46.22   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12711.2  on 9233  degrees of freedom
Residual deviance:  8323.5  on 9232  degrees of freedom
AIC: 8327.5

Number of Fisher Scoring iterations: 5
code R
# Affichage des AIC et deviance
cat("AIC for model with 'identity':", aic_identity, "\n")
AIC for model with 'identity': 11166.46 
code R
cat("AIC for model with 'score':", aic_score, "\n")
AIC for model with 'score': 8327.457 
code R
cat("Deviance for model with 'identity':", deviance_identity, "\n")
Deviance for model with 'identity': 11162.46 
code R
cat("Deviance for model with 'score':", deviance_score, "\n")
Deviance for model with 'score': 8323.457 
code R
library(pROC)

# Calculer la probabilité prédite pour chaque modèle
prob_identity <- predict(model_identity, type = "response")
prob_score <- predict(model_score, type = "response")

# Calculer la courbe ROC et AUC pour chaque modèle
roc_identity <- roc(my_data_frame$Sex, prob_identity)
roc_score <- roc(my_data_frame$Sex, prob_score)

# Afficher les résultats AUC
cat("AUC for model with 'identity':", auc(roc_identity), "\n")
AUC for model with 'identity': 0.7242187 
code R
cat("AUC for model with 'score':", auc(roc_score), "\n")
AUC for model with 'score': 0.8676684 
code R
# Tracer les courbes ROC
plot(roc_identity, col = "blue", main = "Courbes ROC pour 'identity' et 'score'")
plot(roc_score, col = "red", add = TRUE)
legend("bottomright", legend = c("identity", "score"), col = c("blue", "red"), lwd = 2)

code R
# Convertir satisfaction en numérique
my_data_frame <- my_data_frame %>%
  mutate(satisfaction_numeric = case_when(
    satisfaction == "Low" ~ 1,
    satisfaction == "Medium" ~ 2,
    satisfaction == "High" ~ 3,
    TRUE ~ NA_real_  # Gérer les valeurs manquantes si nécessaire
  ))

# Vérifier la conversion
table(my_data_frame$satisfaction_numeric)

   1    2    3 
3503 2605 3117 
code R
# Test de corrélation de Spearman entre satisfaction (ordinale) et distance_quartile (continue)
cor.test(my_data_frame$satisfaction_numeric, my_data_frame$distance_abs_2, method = "spearman")

    Spearman's rank correlation rho

data:  my_data_frame$satisfaction_numeric and my_data_frame$distance_abs_2
S = 1.3324e+11, p-value = 0.07856
alternative hypothesis: true rho is not equal to 0
sample estimates:
        rho 
-0.01831584 
code R
# Convertir 'santé' en numérique (Good = 1, Fair = 2, Poor = 3)
my_data_frame <- my_data_frame %>%
  mutate(santé_numeric = case_when(
    Santé == "Good" ~ 1,
    Santé == "Fair" ~ 2,
    Santé == "Poor" ~ 3,
    TRUE ~ NA_real_
  ))

# Test de corrélation de Spearman entre santé (ordinale) et distance_abs_2
cor.test(my_data_frame$santé_numeric, my_data_frame$distance_abs_2, method = "spearman")

    Spearman's rank correlation rho

data:  my_data_frame$santé_numeric and my_data_frame$distance_abs_2
S = 1.2994e+11, p-value = 0.6901
alternative hypothesis: true rho is not equal to 0
sample estimates:
         rho 
-0.004159822 
code R
# Régression logistique ordinale
library(MASS)
rego_ordinal <- polr(satisfaction ~ distance_abs_2, data = my_data_frame, method = "logistic")
summary(rego_ordinal)
Call:
polr(formula = satisfaction ~ distance_abs_2, data = my_data_frame, 
    method = "logistic")

Coefficients:
                 Value Std. Error t value
distance_abs_2 0.02273    0.03199  0.7105

Intercepts:
           Value    Std. Error t value 
High|Low    -0.6543   0.0340   -19.2524
Low|Medium   0.9511   0.0348    27.3139

Residual Deviance: 20135.53 
AIC: 20141.53 
(9 observations effacées parce que manquantes)

1.6 Sexe ou Genre pour prédire la Satisfaction ?

code R
# Convertir les variables en facteurs si nécessaire
my_data_frame$DIPLOM <- as.factor(my_data_frame$DIPLOM)
my_data_frame$Sex <- as.factor(my_data_frame$Sex)
my_data_frame$CLASSIF <- as.factor(my_data_frame$CLASSIF)
my_data_frame$revenu <- as.factor(my_data_frame$revenu)
my_data_frame$Santé <- as.factor(my_data_frame$Santé)
my_data_frame$satisfaction <- as.factor(my_data_frame$satisfaction)
my_data_frame$SITUA <- as.factor(my_data_frame$SITUA)
my_data_frame$CS2D <- as.factor(my_data_frame$CS2D)

# Créer une fonction pour comparer les modèles avec 'score' et 'Sex' comme prédicteurs
compare_models <- function(variable) {
  # Affichage de la variable actuellement traitée
  cat("Traitement de la variable :", variable, "\n")
  
  # Modèle avec score comme prédicteur
  model_score <- polr(as.formula(paste(variable, "~ score_normalise_2")), data = my_data_frame, method = "logistic")
  
  # Modèle avec sexe comme prédicteur
  model_sex <- polr(as.formula(paste(variable, "~ Sex")), data = my_data_frame, method = "logistic")
  
  # Comparer l'AIC des deux modèles
  aic_score <- AIC(model_score)
  aic_sex <- AIC(model_sex)
  
  # Comparer et retourner le meilleur modèle
  if (aic_score < aic_sex) {
    return(data.frame(variable = variable, best_predictor = "Score", AIC_score = aic_score, AIC_sex = aic_sex))
  } else {
    return(data.frame(variable = variable, best_predictor = "Sex", AIC_score = aic_score, AIC_sex = aic_sex))
  }
}

# Liste des variables d'intérêt à analyser
variables <- c("DIPLOM", "CLASSIF", "SITUA", "satisfaction", "Santé", "revenu", "CS2D")

# Appliquer la fonction pour chaque variable d'intérêt et combiner les résultats
results <- do.call(rbind, lapply(variables, compare_models))
Traitement de la variable : DIPLOM 
Traitement de la variable : CLASSIF 
Traitement de la variable : SITUA 
Traitement de la variable : satisfaction 
Traitement de la variable : Santé 
Traitement de la variable : revenu 
Traitement de la variable : CS2D 
code R
# Afficher les résultats sous forme de tableau
print(results)
      variable best_predictor AIC_score   AIC_sex
1       DIPLOM          Score 41875.234 41877.259
2      CLASSIF            Sex  9541.198  9303.852
3        SITUA            Sex 24323.188 24316.900
4 satisfaction            Sex 20141.821 20141.519
5        Santé            Sex 14570.078 14563.140
6       revenu            Sex 19741.506 19729.137
7         CS2D          Score 29090.088 29114.292

1.7 Genre et Professions

code R
library(readxl)
my_data_frame$CS2D<-as.numeric(my_data_frame$CS2D)
Listes_professions <- read_excel("professions.xls")

mapping_codes_libelles <- Listes_professions[, c("Code", "Libelle")]

my_data_frame <- my_data_frame %>%
  left_join(mapping_codes_libelles, by = c("CS2D" = "Code")) %>%
  mutate(CS2D = Libelle)
# Charger la bibliothèque dplyr pour manipuler les données
library(dplyr)


head(my_data_frame$CS2D)
[1] "Chefs d'entreprise de 10 salariés ou plus"
[2] NA                                         
[3] "Agriculteurs sur petite exploitation"     
[4] NA                                         
[5] NA                                         
[6] NA                                         
code R
summary(my_data_frame$CS2D)
   Length     Class      Mode 
     9234 character character 
code R
# Compter les occurrences des professions
profession_counts <- table(my_data_frame$CS2D)

# Trier les professions par fréquence, du plus grand au plus petit
profession_counts_sorted <- sort(profession_counts, decreasing = TRUE)

# Sélectionner les 10 professions les plus représentées
top_10_professions <- head(profession_counts_sorted, 10)

# Afficher les 10 professions les plus représentées
top_10_professions

                                                Artisans 
                                                     332 
                                   Professions libérales 
                                                     228 
                                Commerçants et assimilés 
                                                     219 
                   Agriculteurs sur moyenne exploitation 
                                                     144 
               Chefs d'entreprise de 10 salariés ou plus 
                                                     111 
                    Agriculteurs sur petite exploitation 
                                                     105 
                  Professeurs, professions scientifiques 
                                                      90 
Professions de l'information, des arts et des spectacles 
                                                      88 
                    Agriculteurs sur grande exploitation 
                                                      79 
                          Cadres de la fonction publique 
                                                      74 
code R
# Filtrer les données pour ne garder que les 10 professions les plus représentées
top_10_data <- my_data_frame[my_data_frame$CS2D %in% names(top_10_professions), ]

# Créer un tableau de proportions
top_10_data <- my_data_frame %>%
  group_by(Sex, CS2D, quartile_distance, score_normalise_2, distance_abs_2, score_scale) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(Sex, CS2D) %>%
  mutate(proportion = count / sum(count))  # Calcul de la proportion


# Séparer les données en deux sous-ensembles (Hommes et Femmes)
df_men <- top_10_data %>% filter(Sex == "Men")
df_women <- top_10_data %>% filter(Sex == "Women")

# Graphique pour les hommes
fig_men <- plot_ly(df_men, 
                   x = ~CS2D, 
                   y = ~proportion, 
                   color = ~quartile_distance, 
                   type = "bar",
                     # Affichage du libellé et proportion
                   textposition = "inside") %>%
  layout(title = "Proportion de score_scale par niveau de diplôme (Hommes)",
         xaxis = list(title = "profession"),
         yaxis = list(title = "Proportion", tickformat = "%"),
         barmode = "stack",  # Empilement des barres
         hovermode = "closest")  # Afficher les informations les plus proches du survol

# Graphique pour les femmes
fig_women <- plot_ly(df_women, 
                     x = ~CS2D, 
                     y = ~proportion, 
                     color = ~quartile_distance, 
                     type = "bar",
                       
                     textposition = "inside") %>%
  layout(title = "Proportion de score_scale par profession (Femmes)",
         xaxis = list(title = "Profession"),
         yaxis = list(title = "Proportion", tickformat = "%"),
         barmode = "stack",  # Empilement des barres
         hovermode = "closest")  # Afficher les informations les plus proches du survol

# Afficher les deux graphiques séparés
fig_men
code R
fig_women
code R
library(dplyr)
library(ggplot2)
library(tidyr)

# Vérifier que CS2D et Sex sont bien des facteurs
top_10_data$CS2D <- as.factor(top_10_data$CS2D)
top_10_data$Sex <- as.factor(top_10_data$Sex)

# Calcul du score moyen par profession
score_summary <- top_10_data %>%
  group_by(CS2D) %>%
  summarise(score_moyen = mean(score_normalise_2, na.rm = TRUE))

# Calcul des proportions Homme/Femme par profession
sex_distribution <- top_10_data %>%
  group_by(CS2D, Sex) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(CS2D) %>%
  mutate(proportion = n / sum(n))

# Fusionner les deux datasets
summary_data <- left_join(score_summary, sex_distribution, by = "CS2D")

# Visualisation combinée
ggplot(summary_data, aes(x = reorder(CS2D, score_moyen))) +
  geom_col(aes(y = score_moyen), fill = "blue", alpha = 0.6) +
  geom_point(aes(y = proportion * max(score_moyen), color = Sex), size = 4) +
  scale_y_continuous(
    sec.axis = sec_axis(~./max(score_summary$score_moyen), name = "Proportion H/F")
  ) +
  labs(title = "Score moyen et proportion H/F par profession",
       x = "Profession (CS2D)", y = "Score normalisé moyen") +
  coord_flip() +
  theme_minimal()

code R
library(dplyr)
library(ggplot2)
library(tidyr)

# Vérifier que CS2D et Sex sont bien des facteurs
top_10_data$CS2D <- as.factor(top_10_data$CS2D)
top_10_data$Sex <- as.factor(top_10_data$Sex)

# Calcul du score moyen par profession
distance_summary <- top_10_data %>%
  group_by(CS2D) %>%
  summarise(distance_moyenne = mean(distance_abs_2, na.rm = TRUE))

# Calcul des proportions Homme/Femme par profession
sex_distribution <- top_10_data %>%
  group_by(CS2D, Sex) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(CS2D) %>%
  mutate(proportion = n / sum(n))

# Fusionner les deux datasets
summary_data <- left_join(distance_summary, sex_distribution, by = "CS2D")

# Visualisation combinée
ggplot(summary_data, aes(x = reorder(CS2D, distance_moyenne))) +
  geom_col(aes(y = distance_moyenne), fill = "blue", alpha = 0.6) +
  geom_point(aes(y = proportion * max(distance_moyenne), color = Sex), size = 4) +
  scale_y_continuous(
    sec.axis = sec_axis(~./max(distance_summary$distance_moyenne), name = "Proportion H/F")
  ) +
  labs(title = "Distance moyenne et proportion H/F par profession",
       x = "Profession (CS2D)", y = "Distance moyenne") +
  coord_flip() +
  theme_minimal()

code R
library(dplyr)
library(ggplot2)

# Vérifier que CS2D et Sex sont bien des facteurs
top_10_data$CS2D <- as.factor(top_10_data$CS2D)
top_10_data$Sex <- as.factor(top_10_data$Sex)

# Calcul de la distance moyenne par profession et par sexe
distance_summary <- top_10_data %>%
  group_by(CS2D, Sex) %>%
  summarise(distance_moyenne = mean(distance_abs_2, na.rm = TRUE), .groups = "drop")

# Visualisation
p<-ggplot(distance_summary, aes(x = reorder(CS2D, distance_moyenne), y = distance_moyenne, fill = Sex)) +
  geom_col(position = "dodge") +
  labs(title = "Distance moyenne (score_distance) par profession et par sexe",
       x = "Profession (CS2D)", y = "Distance moyenne (score_distance)") +
  coord_flip() +
  theme_minimal()

p_plotly<-ggplotly(p)
p_plotly
code R
# Créer un tableau de proportions
top_10_data <- my_data_frame %>%
  group_by(Sex, CS2D, quartile_distance, score_normalise_2, distance_abs_2, score_scale) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(Sex, CS2D) %>%
  mutate(proportion = count / sum(count))  # Calcul de la proportion

library(dplyr)
library(ggplot2)
library(tidyr)

# Vérifier que CS2D, Sex et score_scale sont bien des facteurs
top_10_data$CS2D <- as.factor(top_10_data$CS2D)
top_10_data$Sex <- as.factor(top_10_data$Sex)
top_10_data$score_scale <- as.factor(top_10_data$score_scale)

# Calcul des proportions de score_scale par profession et sexe
score_scale_distribution <- top_10_data %>%
  group_by(CS2D, Sex, score_scale) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(CS2D, Sex) %>%
  mutate(proportion = n / sum(n))

# Visualisation
ggplot(score_scale_distribution, aes(x = score_scale, y = proportion, fill = Sex)) +
  geom_col(position = "dodge") +
  facet_wrap(~ CS2D) +
  labs(title = "Répartition de score_scale par profession et sexe",
       x = "Score Scale", y = "Proportion") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

code R
top_10_data$CS2D <- as.factor(top_10_data$CS2D)
top_10_data$Sex <- as.factor(top_10_data$Sex)
top_10_data$score_scale <- as.factor(top_10_data$score_scale)

# Créer une version numérotée (1-10 + NA) pour affichage dans le graphique
top_10_data$CS2D_numeric <- as.numeric(top_10_data$CS2D)

# Calcul des proportions de score_scale par profession et sexe
score_scale_distribution <- top_10_data %>%
  group_by(CS2D_numeric, CS2D, Sex, score_scale) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(CS2D_numeric, CS2D, Sex) %>%
  mutate(proportion = n / sum(n))

# Graphique ggplot
p <- ggplot(score_scale_distribution, aes(x = score_scale, y = proportion, fill = Sex, 
                                          text = paste("Profession :", CS2D))) +
  geom_col(position = "dodge") +
  facet_wrap(~ CS2D_numeric) +  # Garde les professions sous forme de codes
  labs(title = "Répartition de score_scale par profession et sexe",
       x = "Score Scale", y = "Proportion") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Convertir en ggplotly avec tooltip interactif
ggplotly(p, tooltip = "text")
code R
# Créer une version numérotée (1-10 + NA) pour affichage dans le graphique
top_10_data$CS2D_numeric <- as.numeric(top_10_data$CS2D)

# Calcul des proportions de score_scale par profession et sexe
score_scale_distribution <- top_10_data %>%
  group_by(CS2D_numeric, CS2D, Sex, score_scale) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(CS2D_numeric, CS2D, Sex) %>%
  mutate(proportion = n / sum(n) * 100)  # Conversion en pourcentage

# Graphique ggplot
p <- ggplot(score_scale_distribution, aes(x = score_scale, y = proportion, fill = Sex, 
                                          text = paste("Profession :", CS2D, "<br>",
                                                       "Sexe :", Sex, "<br>",
                                                       "Proportion :", round(proportion, 1), "%"))) +
  geom_col(position = "dodge") +
  facet_wrap(~ CS2D_numeric) +  # Garde les professions sous forme de codes
  labs(title = "Répartition de score_scale par profession et sexe",
       x = "Score Scale", y = "Proportion (%)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Convertir en ggplotly avec tooltip interactif
ggplotly(p, tooltip = "text")